Fit non-gravitational parameters A1/A2/A3 (Yarkovsky) (#351)#356
Open
matthewholman wants to merge 5 commits into
Open
Fit non-gravitational parameters A1/A2/A3 (Yarkovsky) (#351)#356matthewholman wants to merge 5 commits into
matthewholman wants to merge 5 commits into
Conversation
Adds optional fitting of the transverse Marsden non-gravitational parameter A2 (Yarkovsky) alongside the 6D state, the first non-grav parameter layup can solve for. Off by default, so the 6-parameter path is byte-for-byte unchanged. Design (see issue #351): the LM solve is the joint (6+1) system -- state and A2 are correlated and the cross-covariance is the point -- but the 6D state is kept as-is everywhere outside the LM kernel (convert/BK/predict untouched), with A2 carried as an auxiliary parameter. Cartesian engine only for now. ASSIST already supplies d(state)/dA2: a zero-seeded variational particle whose parameter direction dA2=1 is set via assist->particle_params is integrated by the non-grav force routine into the A2 partial. The variational particles are laid out contiguously (var+6 = A2), so the residual/partial loops and the normal-equation assembly simply run to npar = 6 or 7 and read var+i uniformly. - fit_result.cpp: FitResult gains fit_a2, a2, a2_unc (+ bindings). - predict.cpp: add_variational_particles takes n_param extra zero-seed variational particles. - orbit_fit.cpp: compute_residuals[_sequence] enable the non-grav force with the standard 1/r^2 g(r), set particle_params (real A2 + the A2 var direction); compute_single_residuals appends the A2 partial column; compute_dX assembles the npar system; the LM loop steps A2, uses ndof = rows - npar, and reports A2 + its 1-sigma from the joint covariance. run_from_vector_with_initial_guess gains a fit_a2 flag (binding too) and seeds A2 from the initial guess. ASSIST skips the non-grav block when A1=A2=A3=0, so A2 is seeded with a tiny nonzero value to keep its column non-degenerate (the partial is independent of A2's magnitude). Partials use the joint variational Jacobian directly. Validated end-to-end (test_nongrav_a2.py): a 7-parameter fit of a noise-free multi-year arc with a known A2 drives chi-square to ~0 and recovers both the state (~1e-11) and A2 (<1% from a perturbed seed) -- the convergence-to-true-A2 is the check on the partial -- while the 6-parameter fit is left visibly biased. Still TODO (follow-ups): expose through the orbitfit() Python driver (a fit-nongrav flag + A2/A2-unc result columns), a weak-constraint guard, COM/KEP covariance for the joint fit, and A1/A3. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Adds a fit_nongrav flag to orbitfit() (and _orbitfit). When set, after the 6-parameter orbit converges the driver refines it jointly with A2 (seeded from that solution) via run_from_vector_with_initial_guess(fit_a2=True), and appends a2 + a2_unc columns (au/day^2) to the result. Cartesian engine only -- a non-cartesian engine is overridden with a warning, mirroring the radar path. A2 is weakly constrained on short arcs, so if the joint refinement fails to converge the driver keeps the 6-parameter solution and reports A2 as NaN (graceful guard). Schema: a2/a2_unc are added only when fit_nongrav, so the default 6-parameter output dtype is unchanged (_get_result_dtypes / create_empty_result are fit_nongrav-aware). test_nongrav_a2.py adds two driver-level tests: a ground ra/dec arc fit through orbitfit(fit_nongrav=True) recovers A2 to <1% with the new columns present, and the default call produces no a2 columns and still converges (no regression). Follow-ups unchanged from the C++ core commit: weak-constraint guard tuning, COM/KEP joint covariance, A1/A3, and threading an A2 seed from a supplied initial_guess (currently A2 starts near zero and the fit finds it). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
On a short arc the A2 column of the design matrix is nearly collinear with the state directions, so the joint 7-parameter solve is rank-deficient and yields a garbage A2 plus a contaminated, over-confident state. Detect this and refuse it. The signal is the conditioning of the CORRELATION matrix of the normal matrix C (C normalized by its diagonal). Raw rcond(C) is useless here because C mixes AU and AU/day^2 scales; the correlation matrix is unit-invariant and measures the genuine collinearity. When its reciprocal condition number falls below WEAK_NONGRAV_RCOND (1e-8, tunable), orbit_fit returns flag = 6. The orbitfit() driver already falls back to the 6-parameter solution when the A2 refinement does not converge, so a flag-6 result is reported as a clean 6-parameter orbit with a2 = a2_unc = NaN -- no contamination. This is a numerical-degeneracy guard, not a significance cut: a full-rank but loosely constrained A2 (e.g. a few-month arc) is still reported, with an honest (large) a2_unc, so genuine non-detections/upper-limits are preserved. Tests: a ~10-day arc trips flag=6 at the run_from_vector level, and through the driver reports a2=NaN while still recovering the 6-parameter state. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Extends the A2-only non-grav fit to any subset of the Marsden parameters
A1 (radial), A2 (transverse), A3 (normal). Same mechanism, generalized: the
joint LM solves 6 + (number of active params), one zero-seed variational
particle per active param (ASSIST integrates d(state)/dA_i via its
particle_params direction), and the residual/partial loops + normal-equation
assembly already run to npar so they need no per-param branching.
- fit_result.cpp: FitResult carries nongrav_mask + a1/a2/a3 + a1_unc/a2_unc/
a3_unc (replacing the A2-only fields).
- orbit_fit.cpp: bool fit_a2 -> int nongrav_mask (bits 1=A1,2=A2,4=A3) through
orbit_fit / compute_residuals[_sequence] / run_from_vector; an ordered active
list drives particle_params, the per-param column packing, the LM step, and
the reported values/uncertainties (from the joint covariance diagonal). The
weak-constraint guard (flag=6) now applies to any non-grav fit. Fixed the
partial scratch arrays in compute_single_residuals (were sized 7 for A2-only;
now 9 for state + A1/A2/A3, which a 3-param fit overflowed).
- orbitfit.py: fit_nongrav accepts False / True (==A2) / a string or iterable
naming params (e.g. "A1A2A3", ["A1","A3"]); _parse_nongrav -> mask + names;
per-param a{n}/a{n}_unc result columns (only those fitted), so the default
6-parameter schema is unchanged.
Tests (test_nongrav_a2.py): individual A1, A2, A3 recovery; a combined A1+A2+A3
fit recovering all three; the driver "A1A2A3" path; plus the existing guard and
schema-unchanged tests. A1/A3 use larger (test-only) magnitudes since astrometry
constrains them far more weakly than A2. 10/10 pass; 6-parameter regression
(orbit_fit, TNO) unchanged.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
matthewholman
changed the title
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
Adds optional fitting of the non-gravitational Marsden parameters A1 (radial), A2 (transverse / Yarkovsky), A3 (normal) alongside the 6D state — any subset, fit jointly (issue #351). Off by default, so the 6-parameter path is unchanged.
Design
The LM solve is the joint (6+N) system (state + N active non-grav params) — state and the non-gravs are correlated and that cross-covariance is the point — but the 6D state stays as-is everywhere outside the LM kernel (convert / BK / predict untouched), with the params carried alongside. Cartesian engine only.
ASSIST supplies
d(state)/dA_ifor free: one zero-seeded variational particle per active param, whose parameter directiondA_i=1is set viaparticle_params, is integrated by the non-grav force routine. The variational particles are contiguous, so the residual/partial loops and the normal-equation assembly run tonparand readvar+iuniformly — no per-parameter branching.Pieces
FitResultgainsnongrav_mask,a1/a2/a3,a1_unc/a2_unc/a3_unc;run_from_vector_with_initial_guess/orbit_fittake anint nongrav_mask(bits 1=A1, 2=A2, 4=A3); per-param column packing, LM stepping, and reported values + 1-sigma (from the joint covariance). Tiny seed on active params (ASSIST skips the non-grav block when A1=A2=A3=0).flag=6; the driver then falls back to the 6-parameter solution and reports the params as NaN. This is a numerical-degeneracy guard, not a significance cut — a full-rank but loosely constrained param is still reported with an honest (large) uncertainty.orbitfit(fit_nongrav=...)acceptsFalse/True(≡ A2) / a string or iterable naming params ("A2","A1A2A3",["A1","A3"]). For each fitted param it addsa{n}anda{n}_uncresult columns; the default 6-parameter output schema is unchanged.Validation
tests/layup/test_nongrav_a2.py(10 tests): individual A1/A2/A3 recovery, a combined A1+A2+A3 fit recovering all three, the driver multi-param path, the short-arc weak-constraint guard, and the schema-unchanged check. A2 recovers a realistic Apophis-scale value; A1/A3 use larger test-only magnitudes since astrometry constrains them far more weakly. The joint fit drives chi-square to ~0 from a perturbed seed — convergence-to-true-params is the check on the partials. Existing orbit-fit / TNO / real-data / streak / iod_auto suites pass.Motivation
The radar cross-check showed the missing Yarkovsky term dominates long-baseline residuals (Apophis); radar greatly improves non-grav observability, so this and the radar work reinforce each other.
Follow-ups
COM/KEP joint covariance (the params are frame-invariant scalars → block Jacobian), and threading a non-grav seed from a supplied
initial_guess(today they start near zero and the fit finds them).🤖 Generated with Claude Code