From b430d504fa635fc55941a40734ce0c05e186ae2a Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 17:22:29 -0400 Subject: [PATCH] Make the non-grav placeholder negligible and the discarded refinement visible Two small defects in the same path, both of which hide a failed non-grav fit. ASSIST skips the non-gravitational block entirely when A1=A2=A3=0, so a zero seed would leave the parameter's column zero. The fit substitutes a placeholder to avoid that, and the partial does not depend on its magnitude -- but the value chosen, 1e-15 au/day^2, is the order of a real Yarkovsky amplitude. On (6489) Golevka it is 6% of the true A2 and raises the starting chi-square from 12,888 to 138,416 before the first iteration, which is enough to trip the chi-square gate and discard the fit. Use 1e-20, which restores the gravity-only chi-square exactly. When the joint fit does return a non-zero flag the driver keeps the six-parameter solution and reports the non-grav parameters as NaN, which is right, but said so only at DEBUG. The caller sees flag 0, a chi-square belonging to the gravity fit, and NaN parameters, with nothing to explain them. Warn instead, and name the flag the joint fit returned: 6 is a weakly-constrained solution, not a failure to converge, and 2 is a fit rejected on chi-square. Co-Authored-By: Claude Opus 5 --- src/layup/orbitfit.py | 12 +++++++++++- src/lib/orbit_fit/orbit_fit.cpp | 17 +++++++++++++---- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/src/layup/orbitfit.py b/src/layup/orbitfit.py index 32d96b51..adf5267b 100644 --- a/src/layup/orbitfit.py +++ b/src/layup/orbitfit.py @@ -1568,7 +1568,17 @@ def _orbitfit( if res_ng.flag == 0: res = res_ng else: - logger.debug("Non-grav refinement did not converge; reporting non-grav params as NaN.") + # Not silent: the caller asked for a non-grav fit and is getting a + # gravity-only one, with the parameters reported as NaN and a flag + # of 0 that refers to the gravity fit. Name the flag the joint fit + # actually returned -- 6 is a weakly-constrained solution, not a + # failure to converge, and 2 is an acceptable-step fit rejected on + # chi-square. + logger.warning( + "Non-grav refinement returned flag %d; keeping the 6-parameter " + "solution and reporting the non-grav parameters as NaN.", + res_ng.flag, + ) # The non-grav refinement above can replace `res`, so take the check # verdicts from whatever is actually being returned. diff --git a/src/lib/orbit_fit/orbit_fit.cpp b/src/lib/orbit_fit/orbit_fit.cpp index 7a6ca171..d1a0e3f9 100644 --- a/src/lib/orbit_fit/orbit_fit.cpp +++ b/src/lib/orbit_fit/orbit_fit.cpp @@ -93,6 +93,10 @@ namespace orbit_fit // 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; + // Arcseconds per radian (180*3600/pi). Converts astrometric/rate // uncertainties from arcseconds to radians and scales residuals for display. static constexpr double ARCSEC_PER_RAD = 206265.0; @@ -1053,8 +1057,13 @@ namespace orbit_fit // active[] holds the param indices (0=A1,1=A2,2=A3) in column order. ASSIST // skips the non-grav block when A1=A2=A3=0, which would zero the param // columns; seed each active param with a tiny nonzero value so its column - // is non-degenerate (the partial is independent of the param's magnitude, - // and 1e-15 au/day^2 is dynamically negligible). + // is non-degenerate. The partial does not depend on the seed's magnitude, + // so the seed only has to be small enough to be dynamically irrelevant. + // It was 1e-15 au/day^2, which is not: real Yarkovsky amplitudes are of + // that order, so on a long arc with precise data the placeholder is a + // sizeable fraction of the signal. On (6489) Golevka it is 6% of the true + // A2 and raises the starting chi-square from 12,888 to 138,416 before the + // first iteration. 1e-20 restores the gravity-only chi-square exactly. std::vector active; for (int i = 0; i < 3; i++) if (nongrav_mask & (1 << i)) @@ -1079,9 +1088,9 @@ namespace orbit_fit for (int k = 0; k < nactive; k++) { if (a123[active[k]] == 0.0) - a123[active[k]] = 1e-15; + a123[active[k]] = NONGRAV_SEED_PLACEHOLDER; if (per_arc && a123b[active[k]] == 0.0) - a123b[active[k]] = 1e-15; + a123b[active[k]] = NONGRAV_SEED_PLACEHOLDER; } // #419 sequential update: snapshot the prior mean x0 (= the seed state)