diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index 628ef735..66cefbed 100644 --- a/src/lib/orbit_fit/orbit_fit.cpp +++ b/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 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; + } } } diff --git a/tests/layup/test_nongrav_a2.py b/tests/layup/test_nongrav_a2.py index baca9a5f..289402ee 100644 --- a/tests/layup/test_nongrav_a2.py +++ b/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])