Add universal Bernstein-Khushalani engine (bk_native) to do_fit

[matthewholman]

This branch was developed with Claude Code, on top of an initial
design discussion captured in the project notes.

Adds a universal Bernstein-Khushalani (BK) fitter as a new engine
option for do_fit, living entirely inside layup's C++ on top of the
existing variational-integration + residual machinery. No external
dependencies (does not call into liborbfit). The new engine is
always a valid choice -- not distance-dispatched -- because the BK
basis with a fixed bound-orbit energy prior on gdot reduces to
ordinary least-squares for well-arced orbits while smoothly providing
ill-conditioning protection for distant short-arc cases.

What's new

src/lib/orbit_fit/bk_basis.{h,cpp}    BK parameterization primitives (Eigen-only)
src/lib/orbit_fit/bk_fit.cpp          LM driver run_bk_native_fit (reuses
                                      compute_residuals from orbit_fit.cpp)
src/main.cpp                          register bk_basis_bindings + bk_fit_bindings

src/layup/orbitfit.py                 engine='bk_native' dispatch in do_fit
                                      via a new _run_fit(...) helper

tools/bk_engine_sweep.py              full-sweep harness over a diagnostic
                                      dataset; produces a CSV + summary table

tests/layup/test_bk_basis.py          Layer 1: 25 pure-math invariants
tests/layup/test_bk_fit.py            Layer 2: 10 LM-driver smoke + convergence
tests/layup/test_bk_everywhere.py     Layer 3:  6 engine-sweep regressions

How it works (briefly)

BK parameterizes the orbit at epoch by (α, β, γ, α̇, β̇, γ̇) where
γ = 1/r_helio and (α, β) are gnomonic tangent-plane coordinates
of the line-of-sight at a fiducial direction n̂₀. The fit:

1. Pick n̂₀ from the mean of the observations' ρ̂ vectors.
2. Convert the Cartesian seed to BK via cartesian_to_bk.
3. Compute a fixed bound-orbit energy prior on γ̇:
   σ_γ̇² = γ²·(2μγ³ − |α̇·ρ̂_α + β̇·ρ̂_β|²) (this is the exact
   off-fiducial form; the naive γ²(2μγ³ − α̇² − β̇²) is only valid
   at α=β=0 and was caught by an invariant test on first run).
4. LM loop: at each iteration, convert the current BK state to
   Cartesian, call layup's compute_residuals to get tangent-plane
   residuals and the 6-element Cartesian partials, chain-rule through
   the 6×6 dcart_dbk to get BK partials, assemble
   C = BᵀWB + λI + P_prior and grad = BᵀWr + P_prior·p_bk,
   solve, Marquardt-accept/reject, repeat.
5. On convergence, propagate the BK covariance to Cartesian via the
   same 6×6 chain rule and return a standard FitResult.

The Marquardt damping λI and the fixed energy prior P_prior
play complementary roles: λI provides LM's adaptive step control
(large early, decays to zero at convergence), P_prior provides
the persistent physical regularization that survives at convergence
in the line-of-sight direction whenever the data don't constrain it.

Validation on the diagnostic scan

tools/bk_engine_sweep.py runs both engines on every case in a
diagnostic scan (default: 98 cases across 7 populations × 14 arc
lengths, σ_obs = 0.1″, truth state as LM seed). Result:

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 the 79 cases where both engines succeed:

• Drift ratio (BK / Cart): median 0.56, mean 187
• Iteration ratio (BK / Cart): median 0.39, mean 0.52

Headline observations:

• BK never fails when Cartesian succeeds (0/98). In 12 cases
  Cartesian's flag=2 (chi²/dof out of range) while BK converges.
• On the typical case, BK is ~2× closer to truth in ~40% the
  iterations. Mean drift ratio of 187 reflects the distant short-arc
  cases where Cartesian wanders 5-13 AU while BK stays within 0.1 AU.
• The 7 both-failed cases are all 0.04-day (1-hour) arcs with 3
  observations -- a data limitation (0 effective d.o.f.), not an
  algorithm issue.

Standout extreme cases:

• scattered_70AU_arc_014.00d: Cartesian 4.52 AU / 58 iter,
  BK 0.017 AU / 6 iter (≈265× tighter, ≈10× fewer iter).
• sednoid_80AU_arc_007.00d: Cartesian 8.65 AU / 100 iter,
  BK 0.24 AU / 10 iter.

Test coverage (41 tests total, all passing)

• Layer 1 (test_bk_basis.py, 25): Cartesian↔BK round-trip,
  finite-difference vs analytic Jacobian, mixed-partial symmetry,
  fiducial-direction gauge invariance, special-case forms at the
  fiducial, σ_γ̇² agreement with the Cartesian-side parabolic
  boundary. Pure math, no ASSIST needed.
• Layer 2 (test_bk_fit.py, 10): empty-obs and few-obs guards,
  recovery of synthetic mainbelt and TNO orbits from perfect and
  0.1%-perturbed seeds, BK↔Cartesian agreement, engine dispatch via
  _run_fit. Needs the ASSIST ephemeris; skips cleanly if absent.
• Layer 3 (test_bk_everywhere.py, 6): regression tests against
  the diagnostic-scan dataset for well-arced (both engines converge,
  agree, drift < 1% r_helio) and distant short-arc (BK drifts no
  more than Cartesian, uses fewer LM iterations).

Dependencies

Stacks on feat/pybind11-eigen-includes (PR #322 — adds the
pybind11/eigen.h headers required to expose BKState, BKFiducial,
dcart_dbk, etc. to Python). Once #322 lands, this PR collapses to
its 7 new commits.

This branch is intentionally independent of feat/bk-engine-dispatch
(PR #323), which adds engine='auto' for distance-based dispatch
between Cartesian and liborbfit's BK fit. The two PRs can land in
either order; once both are in, engine={'cartesian', 'bk', 'auto', 'bk_native'} are all available and the diagnostic harness can sweep
all four for comparison.

Review Checklist for Source Code Changes

• Does pip install still work?
• Have you written a unit test for any new functions? — 41 tests across test_bk_basis.py, test_bk_fit.py, test_bk_everywhere.py
• Do all the units tests run successfully? — full suite green
• Does Layup run successfully on a test set of input files/databases? — yes, plus full diagnostic-scan sweep (98 cases)
• Have you used black on the files you have updated?

[matthewholman]

Heads up: CI is failing here at the "Install dependencies" step, but the failure is an upstream dependency issue, not anything in this PR.

Root cause: ABI mismatch between the latest assist and rebound wheels on PyPI. The build-isolated environment that pip creates for layup pulls the latest of each, and libassist now references a rebound symbol that the latest librebound no longer exports:

OSError: dlopen(libassist.cpython-312-darwin.so):
  Symbol not found: _reb_simulation_create_from_simulationarchive_with_messages
  Referenced from: libassist.cpython-312-darwin.so
  Expected in:     librebound.cpython-312-darwin.so
CMake Error at CMakeLists.txt:34:
  Could not find assist package. Is assist installed in the build environment?

Reproduces locally on a clean clone today (2026-05-15). Affects every branch, not just this one — main's last successful build was 2026-05-11, and the matrix runs since then have stayed queued or quietly failed. PR #325 passes only because its successful CI ran before the upstream break.

The 41 BK tests in this PR all pass locally against the previously-installed assist 1.2.0 / rebound 4.6.0 combo (i.e., this PR doesn't cause any regression independent of the upstream issue).

Suggested fix (separate from this PR, since it's a repo-wide problem):

[build-system]
requires = [
    "scikit-build-core>=0.10",
    "pybind11",
    "assist>=1.2,<1.3",   # or whichever range is known-compatible
    "rebound>=4.6,<4.7",
]