Report a converged but physically impossible orbit as flag 9 (#493)

src/layup/constants.py

@@ -60,6 +60,32 @@
 FLAG_DEGENERATE_COV = 6  # converged; covariance degenerate, or a variance non-positive
 FLAG_PRIOR_NOT_POSITIVE_DEFINITE = 7  # incremental: prior covariance ill-posed, so a full refit
 FLAG_INCREMENTAL_NO_FULL_OBS = 8  # incremental update with no full observation set to refit from
+FLAG_IMPLAUSIBLE_ORBIT = 9  # converged; hyperbolic excess speed implausibly large
+
+# Hyperbolic excess speed, in km/s, above which a converged fit is reported as
+# physically impossible.
+#
+# A short arc can converge with an excellent reduced chi-square onto a state no
+# real object could occupy, because the arc does not constrain the velocity. The
+# chi-square check cannot catch this: it is anti-correlated with the failure,
+# since the less the object moves across the arc the better the fit.
+#
+# Excess speed is the one statement available without assuming the object is
+# bound. Layup is expected to fit genuine interstellar objects, and those are
+# unbound and fast -- 3I/ATLAS arrives at about 59 km/s -- so boundedness itself
+# is never grounds for rejection, and only a speed far above any plausible
+# arrival speed is evidence of a bad fit rather than an unusual object. The
+# default is about 3.4x the fastest interstellar object yet observed; measured
+# against long-arc truth orbits it rejected no correct fit, and every fit it
+# rejected was wrong.
+#
+# The value means what it says: raise it to accept faster orbits, and set it low
+# to reject anything near-unbound. To switch the check off, set it far above any
+# achievable speed.
+MAX_EXCESS_SPEED_KM_S = 200.0
+
+# km/s expressed in layup's au/day.
+KM_S_IN_AU_PER_DAY = 86400.0 / AU_KM
 
 # How far the fitting pipeline got before it stopped.
 STAGE_NOT_ATTEMPTED = 0

src/layup/orbitfit.py

@@ -40,11 +40,14 @@
     FLAG_CONVERGED,
     FLAG_CSQ_TOO_LARGE,
     FLAG_DEGENERATE_COV,
+    FLAG_IMPLAUSIBLE_ORBIT,
     FLAG_INCREMENTAL_NO_FULL_OBS,
     FLAG_NO_ROOT_CONVERGED,
     FLAG_NO_SOLUTION,
     FLAG_NOT_ATTEMPTED,
     FLAG_PRIOR_NOT_POSITIVE_DEFINITE,
+    KM_S_IN_AU_PER_DAY,
+    MAX_EXCESS_SPEED_KM_S,
     MU_SUN,
     OUTCOME_COLUMNS,
     SPEED_OF_LIGHT,
@@ -55,6 +58,30 @@
     STAGE_NOT_ATTEMPTED,
     STAGE_PRIMARY,
 )
+
+
+def _implausible_excess_speed(state) -> bool:
+    """True when a converged state's hyperbolic excess speed is impossibly large.
+
+    ``state`` is the barycentric equatorial state (au, au/day). Using barycentric
+    coordinates and the Sun's mass rather than the solar system's are both far
+    below the threshold that matters here.
+    """
+    pos = np.asarray(state[:3], dtype=float)
+    vel = np.asarray(state[3:6], dtype=float)
+    if not (np.all(np.isfinite(pos)) and np.all(np.isfinite(vel))):
+        return False  # a NaN state is already reported by the convergence flag
+    r = float(np.linalg.norm(pos))
+    if not r > 0:
+        return False
+    energy = 0.5 * float(np.dot(vel, vel)) - MU_SUN / r
+    if not energy > 0:
+        return False  # bound, or exactly parabolic
+    # bool(), not the numpy scalar the comparison yields: the annotation says
+    # bool, and the value is written to an integer output column.
+    return bool(np.sqrt(2.0 * energy) > MAX_EXCESS_SPEED_KM_S * KM_S_IN_AU_PER_DAY)
+
+
 from layup.convert import convert
 from layup.iod import filter_candidates_by_residual, get_iod, iod_methods, partition_close_approach
 
@@ -1471,10 +1498,18 @@ def _orbitfit(
             else:
                 logger.debug("Non-grav refinement did not converge; reporting non-grav params as NaN.")
 
-        # The non-grav refinement above can replace `res`, so take the gate
+        # The non-grav refinement above can replace `res`, so take the check
         # verdicts from whatever is actually being returned.
         outcome.record(res)
 
+        # Physical plausibility. Applied here, to the result the caller is about
+        # to receive, rather than per candidate: the driver above runs one fit
+        # per initial-orbit candidate and reads a non-zero flag as "this
+        # candidate failed", so a verdict set per candidate is retried away.
+        if res.flag == FLAG_CONVERGED and _implausible_excess_speed(res.state):
+            outcome.failed_physical = True
+            res.flag = FLAG_IMPLAUSIBLE_ORBIT
+
         # Populate our output structured array with the orbit fit results
         success = res.flag == 0
         if not success:

tests/layup/test_implausible_orbit_flag.py

@@ -0,0 +1,129 @@
+"""A fit can converge, satisfy every statistical check, and still describe an
+orbit no object could occupy.
+
+A short arc does not constrain the velocity, so the differential correction can
+settle on an impossibly fast state with an excellent reduced chi-square. The
+chi-square check cannot catch it: the less the object appears to move across the
+arc, the better that fit looks.
+
+The criterion is hyperbolic excess speed, not boundedness. Layup is expected to
+fit genuine interstellar objects, which are unbound and fast -- 3I/ATLAS arrives
+at about 59 km/s -- so being unbound is never itself grounds for rejection.
+"""
+
+import os
+
+import numpy as np
+import pooch
+import pytest
+
+from layup.constants import (
+    FLAG_CONVERGED,
+    KM_S_IN_AU_PER_DAY,
+    MAX_EXCESS_SPEED_KM_S,
+    MU_SUN,
+)
+from layup.orbitfit import _implausible_excess_speed
+
+FASTEST_KNOWN_INTERSTELLAR_KM_S = 59.0  # 3I/ATLAS
+
+
+def _state_with_excess_speed(v_inf_km_s, r_au=1.0):
+    """A state at ``r_au`` whose hyperbolic excess speed is ``v_inf_km_s``.
+
+    From the vis-viva energy, v_inf^2 = v^2 - 2*GM/r, so a speed of
+    sqrt(v_inf^2 + 2*GM/r) at radius r gives exactly the excess speed asked for.
+    """
+    v_inf = v_inf_km_s * KM_S_IN_AU_PER_DAY
+    speed = np.sqrt(v_inf**2 + 2.0 * MU_SUN / r_au)
+    return [r_au, 0.0, 0.0, 0.0, speed, 0.0]
+
+
+# --------------------------------------------------------------------------
+# The criterion itself. Constructed states, so nothing has to move the
+# threshold to exercise both sides of it.
+# --------------------------------------------------------------------------
+
+
+def test_a_bound_orbit_is_never_implausible():
+    """Circular at 1 au: negative energy, so there is no excess speed at all."""
+    speed = np.sqrt(MU_SUN)
+    assert _implausible_excess_speed([1.0, 0.0, 0.0, 0.0, speed, 0.0]) is False
+
+
+def test_just_below_the_threshold_is_accepted():
+    state = _state_with_excess_speed(MAX_EXCESS_SPEED_KM_S * 0.99)
+    assert _implausible_excess_speed(state) is False
+
+
+def test_just_above_the_threshold_is_rejected():
+    state = _state_with_excess_speed(MAX_EXCESS_SPEED_KM_S * 1.01)
+    assert _implausible_excess_speed(state) is True
+
+
+def test_the_threshold_clears_the_fastest_known_interstellar_object():
+    """Rejecting a real interstellar object would reject the bound-to-unbound
+    case layup exists to handle, so the margin is checked rather than assumed."""
+    state = _state_with_excess_speed(FASTEST_KNOWN_INTERSTELLAR_KM_S)
+    assert _implausible_excess_speed(state) is False
+    assert MAX_EXCESS_SPEED_KM_S > 3 * FASTEST_KNOWN_INTERSTELLAR_KM_S
+
+
+def test_a_nan_state_is_left_to_the_convergence_flag():
+    """A diverged fit is already reported as not converged; this check must not
+    turn it into a different kind of failure."""
+    assert _implausible_excess_speed([np.nan] * 6) is False
+
+
+def test_a_state_at_the_origin_is_not_reported():
+    """Radius zero makes the energy undefined rather than large."""
+    assert _implausible_excess_speed([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) is False
+
+
+# --------------------------------------------------------------------------
+# End to end, through a real fit.
+# --------------------------------------------------------------------------
+
+CACHE = str(pooch.os_cache("layup"))
+_EPHEM_AVAILABLE = os.path.exists(os.path.join(CACHE, "linux_p1550p2650.440")) and os.path.exists(
+    os.path.join(CACHE, "sb441-n16.bsp")
+)
+
+requires_ephem = pytest.mark.skipif(
+    not _EPHEM_AVAILABLE,
+    reason=f"ASSIST ephemeris missing at {CACHE}; skipping the end-to-end test.",
+)
+
+
+@requires_ephem
+def test_a_real_interstellar_orbit_is_accepted():
+    """3I/ATLAS is unbound and fast, and layup ships its discovery arc. It must
+    come through with the flag clear and the column clear."""
+    import spiceypy as spice
+    from numpy.lib import recfunctions as rfn
+
+    from layup.orbitfit import _orbitfit
+    from layup.utilities.data_processing_utilities import LayupObservatory
+    from layup.utilities.data_utilities_for_tests import get_test_filepath
+    from layup.utilities.file_io.CSVReader import CSVDataReader
+
+    obs = CSVDataReader(
+        get_test_filepath("3I_ATLAS_ades.csv"), "csv", primary_id_column_name="provID"
+    ).read_rows()
+    helper = LayupObservatory(cache_dir=CACHE)
+    et = np.array([spice.str2et(t) for t in obs["obsTime"]], dtype="<f8")
+    obs = rfn.append_fields(obs, "et", et, usemask=False, asrecarray=True)
+    obs = rfn.merge_arrays(
+        [obs, helper.obscodes_to_barycentric(obs)], flatten=True, asrecarray=True, usemask=False
+    )
+    row = _orbitfit(obs, cache_dir=CACHE, primary_id_column_name="provID", iod="gauss", engine="cartesian")[0]
+
+    state = np.array([row[k] for k in ("x", "y", "z", "xdot", "ydot", "zdot")], dtype=float)
+    energy = 0.5 * np.dot(state[3:], state[3:]) - MU_SUN / np.linalg.norm(state[:3])
+    v_inf = np.sqrt(2 * energy) / KM_S_IN_AU_PER_DAY if energy > 0 else 0.0
+
+    assert v_inf == pytest.approx(
+        FASTEST_KNOWN_INTERSTELLAR_KM_S, abs=5.0
+    ), f"expected an arrival speed near {FASTEST_KNOWN_INTERSTELLAR_KM_S} km/s, got {v_inf:.1f}"
+    assert int(row["flag"]) == FLAG_CONVERGED, f"a real interstellar object was rejected (flag={row['flag']})"
+    assert int(row["failed_physical"]) == 0