Skip to content

Add iod='auto': Gauss with BK-IOD fallback#331

Open
matthewholman wants to merge 14 commits into
mainfrom
feat/iod-auto-fallback
Open

Add iod='auto': Gauss with BK-IOD fallback#331
matthewholman wants to merge 14 commits into
mainfrom
feat/iod-auto-fallback

Conversation

@matthewholman
Copy link
Copy Markdown
Collaborator

This branch was developed with Claude Code, on top of an initial
implementation I had started.

Adds iod="auto" to do_fit: runs Gauss IOD first, falls back to
the universal-BK 5-parameter linear IOD (run_bk_iod) if every
Gauss root fails to seed the LM. iod="gauss" (the default) is
unchanged.

Motivation

The diagnostic-scan sweep in tools/bk_iod_sweep.py (added in
PR #328) shows that Gauss and BK-IOD are complementary.
Across 98 cases, the IOD→LM pipelines succeed on:

IOD converges
Gauss alone 82 / 98
BK-IOD alone 72 / 98
Either 90 / 98

So a single "try Gauss, fall back to BK" path picks up an extra
~8 cases in this dataset. The regimes where BK rescues Gauss are
mostly short-arc distant objects -- exactly where BK-IOD's
parallax-based linear fit handles geometries that Gauss's
three-point method finds degenerate.

What changed

  • do_fit's primary-segment seed loop is refactored from the
    previous "try root 0, elif root 1, elif root 2" chain (which had
    a logic bug -- the second elif could never fire) into a clean
    for-loop over all Gauss roots that breaks on first LM convergence.
  • After the Gauss-roots loop, if iod='auto' and no Gauss root
    succeeded, do_fit calls run_bk_iod on the primary segment
    (using the middle observation's epoch, matching Gauss's
    idx1 convention) and re-runs the LM with that seed.
  • iod='gauss' still raises ValueError on any unknown value;
    iod='auto' is now also accepted.

Tests

tests/layup/test_iod_auto.py (4 tests):

  • test_do_fit_unknown_iod_raises: ValueError on typo'd iod.
  • test_iod_auto_matches_gauss_on_easy_case: on a mainbelt arc
    where Gauss converges, iod='auto' produces a state
    bit-identical to iod='gauss' (the BK fallback never fires).
  • test_iod_auto_recovers_when_gauss_fails: on
    sednoid_80AU_arc_001.00d, iod='gauss' returns flag=3
    (primary interval failed) under the default Cartesian LM,
    while iod='auto' returns flag=0 by falling back to BK-IOD.
  • test_iod_auto_engine_choice_propagates: iod='auto' is
    independent of the engine choice (works with both
    engine='cartesian' and engine='bk_native').

Dependencies

Stacks on feat/bk-iod (PR #328 -- provides run_bk_iod). Once
PR #328 lands, this PR collapses to its 1 new commit.

Review Checklist for Source Code Changes

  • Does pip install still work?
  • Have you written a unit test for any new functions? — tests/layup/test_iod_auto.py
  • Do all the units tests run successfully? — 4 new tests + 56 BK + IOD tests
  • Does Layup run successfully on a test set of input files/databases?
  • Have you used black on the files you have updated?

matthewholman and others added 14 commits May 15, 2026 18:03
… Python

Without these, accessing rho_hat / a_vec / d_vec from Python raises
"Unable to convert function return value to a Python type" because the
binding can't see the Eigen and std::array conversions.

This is a minimal enabler (no behavior change); the A/D-vector
correctness fix on fix/tangent-basis-vectors is independent and lives
on its own branch.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Pure C++/Eigen math layer for the Bernstein-Khushalani parameterization,
with the barycenter as the BK coordinate origin and gnomonic projection
defining the tangent plane at a fiducial direction n0.  No ASSIST,
REBOUND, or pybind11 dependencies, so this translation unit can be
tested in isolation.

New files in src/lib/orbit_fit/:
  bk_basis.h    -- types (BKState, BKFiducial) and function declarations
  bk_basis.cpp  -- implementations: choose_fiducial, bk_to_cartesian,
                   cartesian_to_bk, dcart_dbk (full 6x6 including the
                   bottom-left cross-term block), sigma_gdot_sq

The 6x6 Jacobian dcart_dbk has the block structure expected from the
math:
  [ d r / d (alpha,beta,gamma)     0                          ]
  [ d v / d (alpha,beta,gamma)     d v / d (adot,bdot,gdot)   ]
with the top-left and bottom-right 3x3 blocks identical (both built
from the (1/gamma)-scaled tangent vectors), and the bottom-left block
holding the cross-term contributions through the second derivatives
d^2 rho_hat / d (alpha, beta)^2.

sigma_gdot_sq returns the bound-orbit energy-prior variance,
gamma^2 (2 mu gamma^3 - adot^2 - bdot^2), or +infinity when the
right-hand side is non-positive (tangential rates already exceed
escape).  Returning +infinity yields zero precision in the prior
matrix used by the LM step, which is the desired "no prior" behavior.

orbit_fit.cpp adds a single line to its unity-build chain:
  #include "bk_basis.cpp"

so the new translation unit compiles into the existing _core module.
No pybind11 binding yet -- that comes in a follow-up commit alongside
Python-side unit tests for the math primitives.

The math derivation, design decisions (barycenter origin, fixed
gdot prior, eigendecomp+energy-prior solver, file layout, layered
test plan) live in the project memory file bk_everywhere_design.md.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Adds a bk_basis_bindings(py::module&) entry point that exposes BKState,
BKFiducial, bk_choose_fiducial, bk_to_cartesian, cartesian_to_bk,
dcart_dbk, and sigma_gdot_sq to Python via pybind11 / pybind11/eigen.h,
and wires the binding into main.cpp's _core module alongside the
existing detection_bindings etc.

tests/layup/test_bk_basis.py covers the pure-math invariants
(25 tests, all passing):

  * Round-trip Cartesian <-> BK across mainbelt / NEO / TNO regimes
    (rtol 1e-12).
  * Analytic dcart_dbk vs central-difference, per-element relative
    error < 1e-5 with parameter-scaled epsilon.
  * Mixed-partial symmetry of the second-derivative cross-terms
    appearing in the bottom-left block of dcart_dbk
    (d^2 r / d alpha d beta == d^2 r / d beta d alpha to FD tolerance).
  * Fiducial-direction gauge invariance: two valid n0 choices recover
    the same Cartesian orbit through round-trip.
  * Special-case forms at the fiducial direction alpha = beta = 0:
    position is (1/gamma) n0, top-left and bottom-right Jacobian
    blocks are [(1/gamma) a, (1/gamma) b, -(1/gamma^2) n0] as columns,
    bottom-left block vanishes when rates are zero.
  * sigma_gdot_sq agreement with the Cartesian energy bound at the
    parabolic boundary, and +inf return when tangential rates already
    exceed escape.

The energy-bound test caught a real bug in the first cut of
sigma_gdot_sq: the formula gamma^2 (2 mu gamma^3 - adot^2 - bdot^2)
is only exact at the fiducial direction.  Off-fiducial, the gnomonic
tangent vectors rho_hat_alpha, rho_hat_beta have magnitudes
sqrt((1+beta^2))/s^2 and sqrt((1+alpha^2))/s^2 respectively, and an
inner product -alpha*beta/s^4, so the true tangential-velocity term is

  |adot rho_hat_alpha + bdot rho_hat_beta|^2 =
    [adot^2 (1+beta^2) - 2 adot bdot alpha beta + bdot^2 (1+alpha^2)] / s^4

which reduces to adot^2 + bdot^2 only at alpha = beta = 0.  Fixed in
sigma_gdot_sq (and bk_basis.h documentation) to use the exact form,
which reproduces the parabolic-boundary condition |v|^2 = 2 mu / |r|
to machine precision in the test.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
bk_fit.cpp contains the LM driver that performs an orbit fit in the
universal Bernstein-Khushalani parameterization on top of layup's
existing Cartesian variational machinery.  Included from orbit_fit.cpp
inside the `namespace orbit_fit` block, after the Cartesian helpers
(compute_residuals, create_sequences, get_weight_matrix, converged)
and the Observation/FitResult types are in scope, so no forward
declarations are needed.

The driver structure mirrors run_from_vector_with_initial_guess but
operates in BK basis throughout:

  1. Pick a fiducial direction from the observations' rho_hat vectors
     (mean direction, Gram-Schmidt for the orthonormal a, b).
  2. Convert the Cartesian seed to BK via cartesian_to_bk.
  3. Compute a fixed bound-orbit energy prior precision on gdot from
     the BK seed (1 / sigma_gdot_sq), zero precision otherwise.
  4. LM loop:
       - convert current BK state to Cartesian -> reb_particle
       - call compute_residuals to get tangent-plane residuals and
         Cartesian 6-element partials per observation
       - chain-rule: B_bk = B_cart * dcart_dbk(current BK, fiducial)
       - assemble C = B_bk^T W B_bk + lambda I + P_prior,
                 grad = B_bk^T W r + P_prior * p_bk
       - solve, Marquardt rho-ratio accept/reject, update BK state,
         check convergence (using the existing `converged` predicate).
  5. On convergence:
       - cov_bk = (B^T W B + P_prior)^-1 (Hessian without lambda)
       - cov_cart = J cov_bk J^T (J = dcart_dbk at converged BK state)
       - return FitResult with state = bk_to_cartesian(BK_final) and
         cov flattened from cov_cart.  method = "bk_native".

Initial lambda and Marquardt accept threshold match the Cartesian
fit at orbit_fit.cpp:553.  Early-exit guard: returns a non-success
FitResult (flag = 1) without crashing when detections.size() < 3.

main.cpp gains an orbit_fit::bk_fit_bindings(m) call alongside the
existing orbit_fit bindings, exposing run_bk_native_fit to Python.

tests/layup/test_bk_fit.py covers the Layer 2 smoke tests:
  * binding loads and run_bk_native_fit returns a FitResult
  * empty-obs path returns flag != 0 without crashing
  * <3 obs path triggers the early-exit guard

Layer 2 convergence tests against synthetic observations from a
known orbit (and the Cartesian/BK agreement test on well-arced
mainbelt) are next steps -- they need either the predict-path
output piped back in or the diagnostic/scan dataset wired up to
this branch.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Generate synthetic observations from a known Cartesian state via
layup's predict_sequence path (a fixed barycenter observer, so the
only dynamical content is the orbit itself), then feed those
observations back into both run_bk_native_fit and the existing
Cartesian fit.

Three new test categories on top of the smoke tests:

  * test_bk_native_fit_recovers_known_state: with the truth state as
    the seed, BK converges in essentially one iteration to a fit
    state matching the truth to rtol=1e-6 and chi2 < 1e-12.
    Parameterized over a 3 AU mainbelt 60-day arc and a 40 AU TNO
    300-day arc.

  * test_bk_native_fit_recovers_from_perturbed_seed: with the seed
    perturbed by 0.1% in each component, the LM loop still converges
    to the truth state (rtol=1e-6) -- exercising the chain-rule
    Jacobian + Marquardt damping on a non-trivial number of
    iterations.  Same two orbital regimes.

  * test_bk_and_cartesian_fits_agree: for the well-constrained
    mainbelt case, run_bk_native_fit and run_from_vector_with_initial_guess
    converge to states that agree at rtol=1e-6.  Establishes the
    "no regression on the easy case" baseline.

All seven tests in tests/layup/test_bk_fit.py pass with ASSIST
ephemeris available; skip cleanly without it.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Wires the universal BK fitter (run_bk_native_fit) into layup's Python
do_fit pipeline alongside the existing Cartesian fit, so callers can
choose the engine per call rather than via the C++ entry point.

Changes:

  - src/layup/orbitfit.py
      * Import run_bk_native_fit from layup.routines.
      * Add a module-level _MU_SUN constant (heliocentric GM in
        AU^3 / day^2) used to construct the BK fixed energy prior.
      * Add _run_fit(assist_ephem, initial_guess, obs, engine) helper
        that dispatches to:
          - run_from_vector_with_initial_guess  for engine='cartesian'
          - run_bk_native_fit (with _MU_SUN)    for engine='bk_native'
          - ValueError otherwise
      * do_fit gains an `engine='cartesian'` parameter (default
        preserves the existing behavior).  All five
        run_from_vector_with_initial_guess call sites inside do_fit
        are now routed through _run_fit so the engine choice
        propagates uniformly.

  - tests/layup/test_bk_fit.py
      * test_run_fit_dispatch_cartesian: _run_fit(..., 'cartesian')
        matches direct run_from_vector_with_initial_guess on
        synthetic mainbelt observations.
      * test_run_fit_dispatch_bk_native: _run_fit(..., 'bk_native')
        matches direct run_bk_native_fit(ephem, ig, obs, MU_SUN).
      * test_run_fit_dispatch_unknown_engine_raises: ValueError on
        an unknown engine name.

The 'auto' (distance-dispatched) engine from PR 323 is intentionally
not wired up here; when 323 lands first this branch rebases and
gains both options.  Likewise the 'bk' (liborbfit-backed) engine
from PR 323 is independent of this work.

All 10 tests in test_bk_fit.py and 25 tests in test_bk_basis.py
continue to pass.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
tests/layup/test_bk_everywhere.py drives both engine='cartesian' and
engine='bk_native' against the diagnostic/scan dataset at
~/Dropbox/claude_layup/diagnostic/scan/truth/ -- the same 7-population,
14-arc-length scan that PR 323's auto-dispatch was validated against.
Skips when either the ASSIST ephemeris or the diagnostic scan is
unavailable, so CI is unaffected.

Two test groups:

  * test_engine_sweep_well_arced_cases: on long-arc cases (30-60d
    mainbelt + 60d classical TNO), both engines converge near truth
    (drift < 1% of heliocentric distance) and agree with each other.

  * test_bk_beats_cartesian_on_short_arc_distant: on distant short-arc
    cases (70 AU scattered / 42 AU classical at 10-14 day arcs), BK
    drifts no more than Cartesian from truth AND uses fewer LM
    iterations.  This is the regime BK was designed for, and the
    diagnostic data shows it strongly: on scattered_70AU_arc_014.00d
    BK stays 0.02 AU from truth in 6 iterations while Cartesian
    wanders 4.5 AU over 58 iterations.

The module also exposes a sweep_cases_from_diagnostic() helper for
ad-hoc engine-sweep harness scripts.

All 6 Layer 3 tests pass (in addition to the 25 Layer 1 + 10 Layer 2
tests, for 41 total BK tests on this branch).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
A runnable CLI script that drives both engine='cartesian' and
engine='bk_native' across an entire diagnostic-scan directory, writes
per-case metrics to CSV, and prints a population-level summary
(BK wins / Cartesian wins / per-engine failures / mean iteration
counts, plus median+mean drift and iteration ratios).

Usage:
    python tools/bk_engine_sweep.py --scan-dir <dir> --output <csv>

Defaults discover the project's diagnostic scan at
~/Dropbox/claude_layup/diagnostic/scan/truth and the layup ephemeris
cache at ~/Library/Caches/layup.  Both are overrideable via flags,
so anyone with a compatible truth dataset can reproduce.

Running on the 98-case scan (truth state as LM seed, sigma_arcsec=0.1):

  Population              n   BK win   Cart win   cart fail   bk fail   both fail
  -------------------------------------------------------------------------------
  centaur_15AU           14       13          0           0         0           1
  centaur_25AU           14        9          2           2         0           1
  classical_42AU         14       10          3           0         0           1
  mainbelt_2.5AU         14       10          3           0         0           1
  mainbelt_3.5AU         14       12          1           0         0           1
  scattered_70AU         14        7          2           4         0           1
  sednoid_80AU           14        6          1           6         0           1
  -------------------------------------------------------------------------------
  TOTAL                  98       67         12          12         0           7

  Across 79 cases where both engines succeed:
    drift ratio  (BK / Cart):  median=0.560, mean=187.199
    iter ratio   (BK / Cart):  median=0.386, mean=0.524

Headline: BK never fails when Cartesian succeeds (0 / 98), succeeds in
12 cases where Cartesian flag=2's out, and on the typical case is ~2x
closer to truth in ~40% the iterations.  The mean drift ratio of 187 is
inflated by the 70-80 AU short-arc cases where Cartesian wanders 5-13 AU
while BK stays put.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Pure whitespace -- two blank-line adjustments black wants for PEP 8
spacing.  No test changes; 25 tests still pass.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Layup-side analog of liborbfit's prelim_fit, adapted to barycenter-
origin BK.  Implements a single closed-form 5x5 weighted least-squares
solve over (alpha, beta, gamma, adot, bdot) with gdot pinned to 0.
Uses bk_basis primitives (choose_fiducial, bk_to_cartesian, dcart_dbk,
sigma_gdot_sq).  Per-observation linear model:

    t_i = (jd - (r_obs . n0)/c) - epoch         # light-time corrected
    x_obs = (rho_hat . a) / (rho_hat . n0)      # gnomonic projection
    y_obs = (rho_hat . b) / (rho_hat . n0)
    x_obs = alpha + adot*t_i - gamma*X_i
    y_obs = beta  + bdot*t_i - gamma*Y_i

Returned FitResult.cov is the 6x6 BK covariance (top-left 5x5 from
the fit, (5,5) entry from the bound-orbit energy prior) transformed
to Cartesian via the analytic dcart_dbk Jacobian.

New file: src/lib/orbit_fit/bk_iod.cpp (~180 lines, plus pybind11
binding registered as orbit_fit::bk_iod_bindings).  Inlined into
orbit_fit.cpp inside the namespace block, after bk_fit.cpp.

Smoke tests pass (empty-obs and few-obs early-exit guards).

KNOWN LIMIT (TODO before merge): the model omits the perspective
denominator 1/(1 - gamma*ze).  For TNO-scale gamma*ze ~ 0.024 this
is a ~2-5% multiplicative bias on the recovered gamma -- expected
for a linear IOD (matches liborbfit's prelim_fit) and fine for
SEEDING the LM, but not for stand-alone use.  Three tests in
test_bk_iod.py currently fail because they assume tighter
tolerances than a linear IOD can deliver; see
bk_iod_wip_2026_05_16.md for the design decision pending
between (a) loosen the tolerances + add a "seeds-LM-to-truth"
test, or (b) add a perspective-correction iteration to absorb
the bias.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Rather than trying to make the linear IOD give "final orbit" quality
across all populations (which would require either iteration to absorb
the perspective term or restriction to narrow arc windows), this lands
the IOD as-is and documents what it's actually good for: a seed for
the BK LM fit.

bk_iod.cpp now carries an explicit REGIME OF VALIDITY block based on
an empirical sweep over the full diagnostic-scan dataset.  Summary
(best window over arc length, drift / r_helio):

  population      best arc   best drift
  --------------  --------   -----------
  mainbelt 2.5AU      1 d        1.9%   <- 2d+ catastrophically fails
  mainbelt 3.5AU      1 d        0.9%   <- 2d+ catastrophically fails
  centaur 15 AU     0.5 d        2.5%
  centaur 25 AU     0.5 d        0.4%
  classical 42 AU    10 d        2.1%
  scattered 70 AU     7 d        0.5%
  sednoid 80 AU      10 d        1.0%

The model omits the perspective factor 1/(1 - gamma*ze), so its
inherent accuracy floor scales as ~|gamma * ze| ~ (1/r_helio)*1AU.
That makes it sub-percent for TNOs at sweet-spot arc lengths and
useless for mainbelt at multi-day arcs.  Same regime as liborbfit's
prelim_fit, which also omits the perspective term.  Mainbelt should
use Gauss; BK-IOD is for distant objects.

test_bk_iod.py rewritten to match the documented regime:
  - Smoke tests for empty / few-obs guards (kept).
  - Removed the broken synthetic-observer tests: those used a
    stationary barycenter observer, which gives zero parallax
    baseline and leaves gamma fundamentally unconstrained.
  - test_bk_iod_distant_objects: BK-IOD on TNO sweet-spot cases,
    asserts drift < a few percent of r_helio.
  - test_bk_iod_seeds_lm_to_truth: feed the BK-IOD result into
    run_bk_native_fit and assert the LM converges to within 1% of
    truth.  This is the end-to-end "is BK-IOD useful?" test and
    passes on all four cases tested, including the 60-day classical
    case where the IOD itself is 5% off.

10/10 tests pass.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
A CLI harness that, for each case in a diagnostic scan, runs both
Gauss-IOD (first root only, matching orbitfit.do_fit's typical path)
and BK-IOD, then feeds each into run_bk_native_fit and records the
LM-converged drift from truth + iteration count.  Writes per-case
metrics to CSV and prints a population-level summary table.

Usage:
    python tools/bk_iod_sweep.py --scan-dir <dir> --output <csv>

Defaults discover the project's diagnostic scan at
~/Dropbox/claude_layup/diagnostic/scan/truth and the layup ephemeris
cache at ~/Library/Caches/layup, overrideable via flags.

Results on the 98-case scan (truth obs with sigma_arcsec=0.1, BK LM
applied after each IOD):

  Population              n  G->LM ok  B->LM ok  BK win  G win  BK only  G only  neither
  ---------------------------------------------------------------------------------------
  centaur_15AU           14        12        12      11      1        0       0        2
  centaur_25AU           14        12        12      10      1        1       1        1
  classical_42AU         14        11        13       8      3        2       0        1
  mainbelt_2.5AU         14        12         6       4      1        1       7        1
  mainbelt_3.5AU         14        13         5       4      1        0       8        1
  scattered_70AU         14        10        12       7      2        3       1        1
  sednoid_80AU           14        12        12      10      1        1       1        1
  ---------------------------------------------------------------------------------------
  TOTAL                  98        82        72      54     10        8      18        8

  Across 64 cases where both IOD->LM pipelines converge:
    drift ratio (BK_LM / Gauss_LM):  median=0.257, mean=0.563
    iter ratio  (BK_LM / Gauss_LM):  median=0.677, mean=0.853

Two findings worth highlighting:

  1. The BK-IOD regime-of-validity analysis (in bk_iod.cpp's
     docstring) is validated empirically: BK-IOD seeds the LM to
     tighter fits than Gauss for distant objects (53 wins / 8
     losses across centaur, classical, scattered, sednoid), but
     Gauss wins for mainbelt where BK-IOD fails to seed the LM at
     all in ~half of cases.

  2. BK-IOD and Gauss are COMPLEMENTARY -- BK succeeds on 8 cases
     where Gauss fails; Gauss succeeds on 18 cases where BK fails.
     A pipeline trying one and falling back covers 90 / 98 cases
     vs Gauss alone (82) or BK alone (72).  The 8 remaining "both
     fail" cases are all 0.04-day (1-hour) arcs with only 3 obs --
     a data limitation, not an algorithm gap.

Implementation notes:
  - Only the first Gauss root is tried per case (mirroring
    orbitfit.do_fit, which falls back to roots 1, 2 only when the
    LM on root 0 fails).  Trying all 3 unconditionally made the
    sweep effectively unusable: pathological spurious-root seeds
    burned iter_max=100 with no convergence.
  - Gauss returns its FitResult.flag in an uninitialized state;
    layup's production code (and this script) ignore that field
    and check the LM's flag after.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
When iod='auto', do_fit runs Gauss IOD as usual, tries every Gauss
root in turn as a primary-segment seed, and -- if every Gauss root
fails to converge the LM -- falls back to run_bk_iod() on the
primary segment.  iod='gauss' (the default) is unchanged.

Motivation: the diagnostic-scan sweep in tools/bk_iod_sweep.py
shows Gauss and BK-IOD are *complementary*.  Across 98 cases:
  - Gauss-seeded LM converges on 82
  - BK-IOD-seeded LM converges on 72
  - The union covers 90; only 8 cases (all 0.04-day arcs with
    only 3 obs) defeat both methods.

So a single "try Gauss, fall back to BK" path picks up an extra
~8 cases (in this dataset; the regimes where BK rescues Gauss are
mostly short-arc distant objects -- exactly where BK-IOD's
parallax-based linear fit handles geometries that Gauss's three-
point method finds degenerate).

Implementation:

  - do_fit's primary-segment loop is refactored from the previous
    "try root 0, elif root 1, elif root 2" chain (which had a
    logic bug -- the second elif could never fire, since the first
    elif required len(solns) > 1 AND the elif structure means it
    was already covered by the first if-branch) into a clean for-
    loop over all Gauss roots that breaks on first convergence.
  - After the Gauss-roots loop, if iod='auto' and no root succeeded,
    do_fit calls run_bk_iod on the primary segment with the middle
    observation's epoch (matching the convention Gauss IOD uses)
    and re-runs the LM with that seed.
  - iod='gauss' raises ValueError on any unknown value as before;
    iod='auto' is now also accepted.

tests/layup/test_iod_auto.py (4 tests):

  - test_do_fit_unknown_iod_raises: still ValueError on a typo.
  - test_iod_auto_matches_gauss_on_easy_case: on a mainbelt
    arc where Gauss converges, iod='auto' produces an identical
    state (the BK fallback never fires).
  - test_iod_auto_recovers_when_gauss_fails: on
    sednoid_80AU_arc_001.00d, iod='gauss' returns flag=3
    ('primary interval failed') under Cartesian LM, while
    iod='auto' returns flag=0 by falling back to BK-IOD.
  - test_iod_auto_engine_choice_propagates: iod='auto' is
    independent of the engine choice; works with both
    engine='cartesian' and engine='bk_native'.

All 56 BK + IOD tests pass.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
@mschwamb mschwamb requested a review from Little-Ryugu May 21, 2026 17:29
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant