Reject a non-grav fit only when the covariance cannot support an amplitude

src/lib/orbit_fit/orbit_fit.cpp

@@ -87,12 +87,6 @@ namespace py = pybind11;
 
 namespace orbit_fit
 {
-    // Minimum reciprocal condition number of the normal-matrix CORRELATION matrix
-    // for a non-grav (A2) fit to be considered well constrained (issue #351).
-    // Below this the A2 column is effectively collinear with the state (typical on
-    // short arcs) and the fit is reported as weakly constrained (flag = 6).
-    static constexpr double WEAK_NONGRAV_RCOND = 1e-8;
-
     // Placeholder amplitude (au/day^2) used to keep a non-grav parameter's column
     // non-degenerate when its seed is zero; see the seeding loop in orbit_fit.
     static constexpr double NONGRAV_SEED_PLACEHOLDER = 1e-20;
@@ -1323,31 +1317,37 @@ namespace orbit_fit
 
         cov = C.inverse();
 
-        // Weak-constraint guard for the non-grav fit (issue #351). On a short arc
-        // a non-grav column becomes nearly collinear with the state directions, so
-        // the joint solve yields garbage params and a contaminated state. Detect
-        // this via the conditioning of the CORRELATION matrix of the normal matrix
-        // C -- the raw rcond(C) is useless here because C mixes AU and AU/day^2
-        // scales, but the correlation matrix (C normalized by its diagonal) is
-        // unit-invariant and measures genuine collinearity. If it is near-singular,
-        // mark the fit weakly constrained (flag = 6) so the driver falls back to
-        // the 6-parameter solution and reports the non-grav params as NaN.
+        // Degeneracy guard for the non-grav fit (issue #351). Report the fit as
+        // weakly constrained (flag = 6), so the driver falls back to the
+        // six-parameter solution, only when the covariance cannot support an
+        // amplitude at all: a non-positive or non-finite variance means the
+        // inverse has lost that column.
+        //
+        // This used to reject on the conditioning of the whole npar x npar
+        // correlation matrix, which is not a statement about the non-grav column.
+        // A six-parameter orbit fit is itself strongly correlated, so the test was
+        // dominated by the orbit underneath: on (6489) Golevka the state block
+        // alone has rcond 2.6e-10, twenty-five times below the old 1e-8 threshold,
+        // while A2 correlates with the state at only 0.91 to 0.93. Radar tightens
+        // the state and lowers that rcond further, so adding seven radar
+        // observations to 2062 Aten's 980 optical ones turned a reported
+        // A2 = -14.96 +/- 1.70 into NaN. Better data made the fit be rejected.
+        //
+        // The concern the old test was standing in for -- a contaminated,
+        // over-confident amplitude -- was a symptom of the step being taken from
+        // the squared-condition normal matrix, and does not survive that fix. A
+        // ten-day arc now returns A2 with an uncertainty five million times its
+        // own value, which tells the caller exactly what it needs to know, and
+        // unlike a significance test it does not suppress a well-constrained null.
         if (nactive > 0 && flag == 0)
         {
-            Eigen::VectorXd diag = C.diagonal();
-            if ((diag.array() > 0.0).all())
+            for (int j = 6; j < npar; j++)
             {
-                Eigen::VectorXd s = diag.array().rsqrt(); // 1/sqrt(C_ii)
-                Eigen::MatrixXd R = s.asDiagonal() * C * s.asDiagonal();
-                Eigen::JacobiSVD<Eigen::MatrixXd> svd(R);
-                const Eigen::VectorXd &sv = svd.singularValues();
-                double rcond = sv(sv.size() - 1) / sv(0);
-                if (!(rcond >= WEAK_NONGRAV_RCOND)) // also catches NaN
+                if (!(cov(j, j) > 0.0) || !std::isfinite(cov(j, j)))
+                {
                     flag = 6;
-            }
-            else
-            {
-                flag = 6; // non-positive variance -> degenerate
+                    break;
+                }
             }
         }
 

tests/layup/test_nongrav_a2.py

@@ -193,14 +193,45 @@ def test_six_param_fit_unaffected_and_biased_by_a2():
     assert fit6.csq > 1e3 * fit7.csq
 
 
-def test_a2_weak_constraint_guard_on_short_arc():
-    """On a short arc the A2 column is nearly collinear with the state, so the
-    joint fit is rank-deficient. The fitter must flag this (flag=6) rather than
-    returning a contaminated, over-confident solution."""
+def test_a2_on_short_arc_is_reported_as_undetermined():
+    """A ten-day arc cannot separate A2 from the state. The fit must say so in the
+    uncertainty rather than by refusing: the amplitude comes back with a formal
+    error orders of magnitude larger than itself, and larger than any physical
+    Yarkovsky amplitude, which is what tells the caller it is not a detection.
+
+    This is a stronger statement than the flag=6 this used to assert. Refusing
+    only says the fitter declined; an honest uncertainty says how badly the arc
+    constrains A2, and it does not suppress a well-constrained null."""
     obs = _build_arc(arc_days=10.0, n=10)  # ~10 days: A2 not separable from state
     ephem = get_ephem(CACHE)
     fit = run_from_vector_with_initial_guess(ephem, _seed(_STATE), obs, 100, _BIT["A2"])
-    assert fit.flag == 6, f"expected weak-constraint flag 6, got {fit.flag}"
+    assert fit.flag == 0, f"expected a converged fit, got flag {fit.flag}"
+    assert np.isfinite(fit.a2_unc) and fit.a2_unc > 0
+    # Yarkovsky amplitudes for near-Earth asteroids run to ~1e-13 au/day^2; the
+    # formal error here is many orders above that, and above |A2| itself.
+    assert fit.a2_unc > 1e-11, f"expected an uninformative sigma, got {fit.a2_unc:.3e}"
+    assert fit.a2_unc > 1e3 * abs(fit.a2)
+
+
+def test_a2_uncertainty_tracks_arc_length():
+    """The counterpart to the test above: the reported uncertainty is what carries
+    the information, and it improves by orders of magnitude with arc length. Over
+    four years the amplitude is recovered to well under a percent, where over ten
+    days it is not recovered at all -- so the short-arc result is a statement
+    about the arc, not a blanket refusal to report."""
+    ephem = get_ephem(CACHE)
+    short = run_from_vector_with_initial_guess(
+        ephem, _seed(_STATE), _build_arc(arc_days=10.0, n=10), 100, _BIT["A2"]
+    )
+    long = run_from_vector_with_initial_guess(
+        ephem, _seed(_STATE), _build_arc(arc_days=4 * 365.0, n=36), 100, _BIT["A2"]
+    )
+    assert short.flag == 0 and long.flag == 0
+    # The long arc recovers the injected amplitude; the short one does not.
+    assert abs(long.a2 - _TRUE_A2) < 0.01 * abs(_TRUE_A2)
+    assert abs(short.a2 - _TRUE_A2) > abs(_TRUE_A2)
+    # And says so: five orders of magnitude between the two formal errors.
+    assert long.a2_unc < 1e-5 * short.a2_unc
 
 
 def _build_arc_array(arc_days=4 * 365.0, n=36, a123=(0.0, _TRUE_A2, 0.0)):
@@ -337,15 +368,17 @@ def test_orbitfit_default_schema_unchanged():
     assert fit[0]["flag"] == 0
 
 
-def test_orbitfit_driver_short_arc_reports_a2_nan():
-    """Through the driver, a short arc that cannot constrain A2 falls back to the
-    6-parameter solution: flag is success, the state is still recovered, and a2 is
-    reported as NaN (not a contaminated value)."""
+def test_orbitfit_driver_short_arc_reports_an_uninformative_a2():
+    """Through the driver, a short arc that cannot constrain A2 still reports one,
+    with an uncertainty that shows it is not a detection. The state is unaffected:
+    the point of the earlier NaN was to protect the caller from a contaminated
+    orbit, and the orbit is still clean."""
     data, guess = _build_arc_array(arc_days=10.0, n=10)
     fit = orbitfit(data, cache_dir=CACHE, initial_guess=guess, fit_nongrav=True)
     row = fit[0]
-    assert row["flag"] == 0  # 6-parameter fallback succeeded
-    assert np.isnan(row["a2"]) and np.isnan(row["a2_unc"])
+    assert row["flag"] == 0
+    assert np.isfinite(row["a2"]) and np.isfinite(row["a2_unc"])
+    assert row["a2_unc"] > 1e3 * abs(row["a2"])
     # The fallback state is the clean 6-parameter orbit, still near truth.
     state = np.array([row["x"], row["y"], row["z"], row["xdot"], row["ydot"], row["zdot"]])
     pos_rel = np.linalg.norm(state[:3] - _STATE[:3]) / np.linalg.norm(_STATE[:3])