hyperjet 1.14.0

//! Gain-ratio Levenberg-Marquardt driver.
//!
//! Minimizes \\( \Phi(\mathbf{x}) = \sum_i r_i(\mathbf{x})^2 \\)
//! (optionally plus a Gaussian prior penalty
//! \\( (\mathbf{x}-\mathbf{x}_0)^\top P_0^{-1} (\mathbf{x}-\mathbf{x}_0) \\))
//! with a trust-region-style accept/reject loop:
//!
//! 1. One full system/Jacobian evaluation per **outer** iteration, at
//!    accepted points only.
//! 2. An **inner** trial loop solves the damped system
//!    \\( (A + \mu D^2)\,\mathbf{h} = \mathbf{g} \\) and evaluates a
//!    cheap trial cost at \\( \mathbf{x} + \mathbf{h} \\). A trial that
//!    fails the gain-ratio acceptance test is **rejected**: the iterate
//!    stays put, \\( \mu \\) is raised, and the system is re-solved with
//!    the *same* Jacobian.
//! 3. The damping parameter follows Nielsen's rule, driven by the gain
//!    ratio \\( \rho = (\Phi - \Phi_t) / \mathrm{pred} \\) — the ratio
//!    of actual to model-predicted cost reduction.
//!
//! State (the iterate, the held system, problem-side commits) moves
//! **only on acceptance**, so the returned solution — cost, covariance,
//! and problem diagnostics — always corresponds to the returned
//! \\( \mathbf{x} \\), and the cost is monotone non-increasing across
//! accepted iterates.
//!
//! The solver is deterministic by construction: fixed-order scalar
//! `f64` arithmetic only (`+`, `-`, `*`, `/`, `sqrt`, comparisons), no
//! libm transcendentals, no randomness, no threading. Identical inputs
//! produce identical bits on every IEEE 754 platform.
//!
//! Two entry points share one driver core:
//!
//! - [`solve`] — residual-level problems ([`ResidualProblem`]): the
//!   driver assembles the normal equations from residual/Jacobian rows
//!   and owns the optional prior.
//! - [`solve_system`] — normal-equations-level problems
//!   ([`SystemProblem`]): the problem hands the driver an assembled
//!   \\( (\Phi, A, \mathbf{g}) \\) triple per evaluation (e.g. a
//!   Schur-reduced system); the problem owns its complete objective,
//!   priors included.
//!
//! # References
//!
//! - Madsen, K., Nielsen, H.B. & Tingleff, O. (2004), *Methods for
//!   Non-Linear Least Squares Problems*, 2nd ed., IMM/DTU — Algorithm
//!   3.16, eqs. 2.18, 2.21, 3.14-3.15.
//! - Nielsen, H.B. (1999), *Damping Parameter in Marquardt's Method*,
//!   IMM-REP-1999-05, DTU.
//! - Moré, J.J. (1978), "The Levenberg-Marquardt Algorithm:
//!   Implementation and Theory", *Numerical Analysis*, LNM 630 —
//!   running-max column-norm scaling, best-retained-point invariant.
//! - Moré, Garbow & Hillstrom (1980), *User Guide for MINPACK-1*,
//!   ANL-80-74 — one-Jacobian-per-outer-iteration structure,
//!   ftol/xtol/gtol termination, acceptance threshold 1e-4.
//! - Marquardt, D.W. (1963), SIAM J. Appl. Math. 11(2);
//!   Levenberg, K. (1944), Q. Appl. Math. 2.
//!
//! # Example
//!
//! ```
//! use hyperjet::optimization::lm::{solve_nlls, LMConfig, NLLSEvaluation};
//!
//! #[derive(Debug)]
//! struct NoError;
//! impl std::fmt::Display for NoError {
//!     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
//!         write!(f, "no error")
//!     }
//! }
//! impl std::error::Error for NoError {}
//!
//! // Fit y = a*x + b to (0,1), (1,3), (2,5), (3,7)  =>  a = 2, b = 1.
//! let xs = [0.0_f64, 1.0, 2.0, 3.0];
//! let ys = [1.0_f64, 3.0, 5.0, 7.0];
//! let residual = |p: &[f64; 2], x: f64, y: f64| p[0] * x + p[1] - y;
//!
//! let solution = solve_nlls(
//!     |p: &[f64; 2]| {
//!         let residuals: Vec<f64> =
//!             xs.iter().zip(&ys).map(|(&x, &y)| residual(p, x, y)).collect();
//!         let jacobian: Vec<[f64; 2]> = xs.iter().map(|&x| [x, 1.0]).collect();
//!         let cost = residuals.iter().map(|r| r * r).sum();
//!         Ok::<_, NoError>(NLLSEvaluation { residuals, jacobian, cost })
//!     },
//!     |p: &[f64; 2]| {
//!         Ok::<_, NoError>(
//!             xs.iter()
//!                 .zip(&ys)
//!                 .map(|(&x, &y)| {
//!                     let r = residual(p, x, y);
//!                     r * r
//!                 })
//!                 .sum(),
//!         )
//!     },
//!     [0.0; 2],
//!     &LMConfig::default(),
//!     None,
//! )
//! .unwrap();
//!
//! assert!(solution.converged);
//! assert!((solution.x[0] - 2.0).abs() < 1e-8);
//! assert!((solution.x[1] - 1.0).abs() < 1e-8);
//! ```

// Fixed-order indexed loops are deliberate throughout this module:
// the accumulation ORDER is part of the bit-determinism contract, and
// the index-based form keeps that order explicit.
#![allow(clippy::needless_range_loop)]

use crate::linalg::generic::{mat_cholesky, mat_inv, mat_symmetric_eigen, mat_symmetrize};

/// Result of evaluating the residual function at a point.
///
/// The caller pre-weights residuals and Jacobian by the observation
/// weight Cholesky factor: \\(\mathbf{r}_w = L\mathbf{r}\\),
/// \\(J_w = LJ\\) where \\(W = LL^T\\). The solver then accumulates
/// \\(J^T J\\) and \\(J^T r\\) without an explicit weight matrix.
#[derive(Clone, Debug)]
pub struct NLLSEvaluation<const N: usize> {
    /// Pre-weighted residual values.
    pub residuals: Vec<f64>,
    /// Pre-weighted Jacobian rows: \\(\partial r_i / \partial x_j\\).
    /// Must have the same length as `residuals`.
    pub jacobian: Vec<[f64; N]>,
    /// Total cost \\(\sum r_i^2\\). Used for Levenberg-Marquardt adaptation.
    pub cost: f64,
}

/// Optional Bayesian prior on the parameters.
///
/// Augments the normal equations: \\(N \leftarrow N + P_0^{-1}\\),
/// \\(d \leftarrow d - P_0^{-1}(\mathbf{x} - \mathbf{x}_0)\\).
#[derive(Clone, Debug)]
pub struct NLLSPrior<const N: usize> {
    /// Prior mean.
    pub mean: [f64; N],
    /// Inverse of prior covariance matrix.
    pub covariance_inv: [[f64; N]; N],
}

// ── Tuning constants ────────────────────────────────────────────────

/// Reseed floor applied to \\(\mu\\) on rejection:
/// \\(\mu \leftarrow \max(\mu, \mu_{\text{seed}}) \cdot \nu\\).
/// Guarantees no value of \\(\mu\\) (including 0) is an absorbing
/// state under the multiplicative escalation.
const MU_SEED: f64 = 1e-12;

/// Relative floor on the Moré scaling diagonal:
/// \\(d_j \ge 10^{-12} \max_k d_k\\). Relative (not absolute) so the
/// floor is invariant under a global rescaling of the parameters —
/// the parameter vector mixes units spanning many orders of magnitude.
const D_FLOOR_REL: f64 = 1e-12;

// ── Configuration ───────────────────────────────────────────────────

/// Solver configuration.
///
/// Construct via [`Default`] and mutate fields; the struct is
/// `#[non_exhaustive]`, so struct-literal construction is reserved to
/// this crate and new knobs can be added without breaking consumers.
#[derive(Clone, Debug)]
#[non_exhaustive]
pub struct LMConfig {
    /// Maximum **outer** iterations (each holds one system/Jacobian
    /// snapshot). Must be ≥ 1.
    ///
    /// This bounds ITERATIONS, not evaluations. Full evaluations
    /// (residuals + Jacobian) =
    ///
    /// \\[
    ///   1_{\text{initial}} \;+\; \text{accepted steps} \;+\;
    ///   \text{rolled-back acceptances} \;+\;
    ///   \text{model refreshes}
    /// \\]
    ///
    /// The last term is the
    /// [`refresh_model`](CostProblem::refresh_model) hook: every `true`
    /// declaration costs one additional full evaluation, reported as
    /// [`LMSolution::n_model_refreshes`]. It is offered at most once
    /// per outer iteration, so that term is bounded by
    /// [`LMSolution::iterations`] — but it is NOT charged against this
    /// budget, and a solve that refreshes on every iteration performs
    /// up to `max_iterations` full evaluations more than the first
    /// three terms imply. Problems that never refresh (the default
    /// hook returns `false`) spend none.
    pub max_iterations: usize,
    /// Trial budget per outer iteration (rejections re-use the same
    /// system). Must be ≥ 1.
    pub max_inner_trials: usize,
    /// Initial damping scale:
    /// \\(\mu_0 = \tau \cdot \max_j (A_{jj} / d_j^2)\\)
    /// (Nielsen eq. 3.14, GSL's scaled variant). Must be finite
    /// and > 0. 1e-3 is the general default; use ~1e-6 for warm
    /// restarts near a converged fit and ~1.0 for poor seeds —
    /// although the accept/reject loop makes the solver insensitive
    /// to this choice (see the τ-insensitivity test).
    pub tau: f64,
    /// Damping cap: exceeding it while rejecting terminates with
    /// [`TerminationReason::DampingExhausted`]. Must be finite and > 0.
    pub mu_max: f64,
    /// Acceptance threshold \\(\eta\\): accept iff
    /// \\(\Phi - \Phi_t > \eta \cdot \mathrm{pred}\\) (MINPACK's 1e-4).
    /// Must be finite and ≥ 0.
    pub min_relative_decrease: f64,
    /// Scaled-gradient tolerance (MINPACK cosine form): converge when
    /// \\(|g_j| \le \texttt{gtol} \cdot d_j \sqrt{\Phi}\\) for all
    /// \\(j\\). 0 disables. Must be finite and ≥ 0.
    pub gtol: f64,
    /// Quadratic-form step tolerance (Milani \\(\lVert\cdot\rVert_C\\)
    /// convention): converge when
    /// \\(\mathbf{h}^\top A \mathbf{h} \le \texttt{qtol}\\) under the
    /// undamped, prior-augmented \\(A\\) that produced the step.
    /// 0 disables (the default — opt in for differential-correction
    /// style consumers). Must be finite and ≥ 0.
    pub qtol: f64,
    /// Relative step tolerance (Nielsen eq. 3.15b): converge when
    /// \\(\lVert\mathbf{h}\rVert_D \le \texttt{xtol}\,(\lVert\mathbf{x}\rVert_D + \texttt{xtol})\\).
    /// Degenerate in delta-coordinate formulations where
    /// \\(\lVert\mathbf{x}\rVert \approx 0\\) — use `qtol` there.
    /// 0 disables. Must be finite and ≥ 0.
    pub xtol: f64,
    /// Relative cost-reduction tolerance (MINPACK `ftol` analogue):
    /// converge when both the actual and the predicted reduction of the
    /// last accepted step are ≤ `ftol`·Φ and the gain ratio of that
    /// step is ≤ 2. 0 disables. Must be finite and ≥ 0.
    pub ftol: f64,
    /// Consecutive invalid (`Err` or non-finite) trial costs before the
    /// solve fails with [`LMError::PersistentInvalidTrials`]. Must
    /// be ≥ 1.
    pub max_consecutive_invalid: usize,
    /// Geodesic acceleration (Transtrum & Sethna 2012; GSL
    /// `multifit_nlinear`): augment each trial step with the
    /// second-order correction \\(\mathbf{h} = \mathbf{v} +
    /// \tfrac{1}{2}\mathbf{a}\\), where \\((A + \mu D^2)\,
    /// \mathbf{a} = -J^\top \mathbf{r}''_{vv}\\) and
    /// \\(\mathbf{r}''_{vv}\\) is the problem-supplied directional
    /// second derivative of the residuals along the velocity step
    /// (see [`CostProblem::second_directional_derivative`]).
    /// Dramatically reduces iteration counts on curved-valley
    /// ("sloppy") cost surfaces. Default `false`; also inert while
    /// the problem returns `None` from the hook or on the
    /// normal-equations path (no Jacobian rows).
    pub geodesic_acceleration: bool,
    /// Correct the predicted reduction with the EXACT second-order term
    /// \\(h^\\top S h\\), \\(S = \\sum_i r_i \\nabla^2 r_i\\)
    /// (default: `false`).
    ///
    /// Gauss-Newton approximates \\(\\nabla^2\\Phi = J^\\top J + S\\) by
    /// \\(J^\\top J\\). Dropping \\(S\\) is harmless only when the
    /// residuals become small at the solution. On a LARGE-RESIDUAL fit it is
    /// not: near the floor \\(J^\\top r \\to 0\\), so the linearized
    /// problem reports itself stationary and the predicted reduction
    /// collapses, while the true objective still has descent available. The
    /// gain ratio then runs to hundreds and the consistency guard refuses
    /// convergence — measured on a joint optical + radar orbit fit at
    /// \\(\\chi^2/\\nu \\approx 13\\), where \\(\\rho\\) went 0.99, 65,
    /// 297, 1061 across four iterations.
    ///
    /// NL2SOL (Dennis, Gay & Welsch 1981) carries a SECANT approximation to
    /// \\(S\\). This does not approximate it. \\(h^\\top S h =
    /// \\sum_i r_i (h^\\top \\nabla^2 r_i h)\\) is a residual-weighted
    /// second directional derivative along \\(h\\) — exactly what
    /// [`ResidualProblem::second_directional_derivative`] returns — so the
    /// correction costs ONE second-order evaluation per trial and \\(S\\) is
    /// never formed.
    ///
    /// This corrects the MODEL, not the step: the step still comes from the
    /// Gauss-Newton system. It makes the gain ratio measure what it is
    /// defined to measure, which is what a large-residual fit needs before
    /// any of LM's acceptance and convergence logic is meaningful.
    ///
    /// Requires the residual path and a problem that implements the
    /// second-derivative hook; declines silently to the Gauss-Newton model
    /// when either is absent, since a missing second derivative is a
    /// capability gap and not an error.
    pub exact_second_order_model: bool,
    /// Compute the step by factoring the least-squares system directly
    /// instead of forming and factoring \(A + \mu D^2\) (default:
    /// `false`).
    ///
    /// Forming the normal equations squares the condition number,
    /// \(\kappa(A^\top A) = \kappa(A)^2\). When observation rows carry
    /// weights spanning orders of magnitude — a joint orbit fit mixing
    /// radar delay at microsecond \(\sigma\) with optical astrometry at
    /// arcsecond \(\sigma\) differs by five orders in
    /// \(\lvert J\rvert/\sigma\) — that squaring can put the system past
    /// f64 range while the underlying problem is perfectly well posed
    /// (Van Loan 1985). Damping then also enters as appended
    /// \(\sqrt{\mu}\) rows rather than as \(\mu D^2\) on a squared
    /// matrix, so large \(\mu\) cannot overflow the step.
    ///
    /// Only the STEP changes. The gain ratio, gradient test, predicted
    /// reduction and reported covariance stay defined on \(A\) and
    /// \(g\), so acceptance behaviour is unchanged on any problem where
    /// both paths can compute a step at all.
    ///
    /// Requires the residual path (rows exist) and is inert on the
    /// system path. Geodesic acceleration is unavailable here — it needs
    /// the Cholesky factor of the damped normal matrix — and declines
    /// rather than silently mixing solvers.
    pub square_root_solve: bool,
    /// Acceleration acceptance guard (GSL `avmax`): the accelerated
    /// step is used only when \\(\lVert\mathbf{a}\rVert_D /
    /// \lVert\mathbf{v}\rVert_D \le \texttt{avmax}\\); beyond it
    /// the truncated expansion is untrustworthy and the trial is
    /// rejected through the normal μ escalation. Must be finite and
    /// > 0 (default 0.75, GSL's default).
    pub avmax: f64,
}

impl Default for LMConfig {
    fn default() -> Self {
        Self {
            max_iterations: 100,
            max_inner_trials: 30,
            tau: 1e-3,
            mu_max: 1e32,
            min_relative_decrease: 1e-4,
            gtol: 1e-8,
            qtol: 0.0,
            xtol: 1e-8,
            ftol: 1.49e-8,
            max_consecutive_invalid: 5,
            geodesic_acceleration: false,
            avmax: 0.75,
            // Off until measured per problem class: it changes the step
            // on any system where the two paths disagree numerically.
            square_root_solve: false,
            // Costs one second-order evaluation per trial; off until a
            // problem class is measured to need it.
            exact_second_order_model: false,
        }
    }
}

// ── Problem traits ──────────────────────────────────────────────────

/// Trial-cost and step-lifecycle interface shared by both problem
/// levels.
pub trait CostProblem<const N: usize> {
    /// Domain error type, propagated through [`LMError`] without
    /// flattening.
    type Error: std::error::Error + 'static;

    /// Objective value at a **trial** point.
    ///
    /// Contract: the same objective, over the same fixed observation
    /// set, as the full evaluation — and no commits to persistent
    /// state (stage into pending slots; commit in
    /// [`on_step_accepted`](Self::on_step_accepted)).
    ///
    /// `Err` or a non-finite value rejects the trial through the
    /// normal damping escalation (bounded by
    /// [`LMConfig::max_consecutive_invalid`]); it is not fatal.
    fn evaluate_cost(&mut self, x: &[f64; N]) -> Result<f64, Self::Error>;

    /// Optional problem-driven clamp on a proposed step.
    ///
    /// The driver computes the predicted reduction from the **actual**
    /// (possibly clamped) step, so the clamp stays visible to the gain
    /// ratio; a clamped step can never declare convergence by itself.
    fn constrain_step(&mut self, _x: &[f64; N], _delta: &mut [f64; N]) {}

    /// Commit hook: persistent state (caches, baselines, diagnostics)
    /// moves ONLY here. Called once for the validated initial
    /// evaluation at \\(\mathbf{x}_0\\) and then after each full
    /// evaluation at an accepted point has succeeded, so committed
    /// state always corresponds to a retained iterate — including on
    /// zero-acceptance exits. Also called after the re-assembly a
    /// [`refresh_model`](Self::refresh_model) declaration triggers,
    /// which lands on the CURRENT accepted point rather than a new one.
    fn on_step_accepted(&mut self, _x: &[f64; N]) {}

    /// Discard hook: drop pending (staged) state from a rejected trial
    /// evaluation. Diagnostics and counters only — nothing committed
    /// needs rolling back.
    fn on_step_rejected(&mut self, _x_trial: &[f64; N]) {}

    /// Model-refresh hook for problems that evaluate an **inexact**
    /// model — a surrogate anchored at a reference point, a cached
    /// linearization, a reduced-fidelity dynamics. Called once per
    /// outer iteration at the current accepted point, after that
    /// iteration's convergence tests and before its first trial.
    ///
    /// This is the ONLY point at which such a problem may change the
    /// model its evaluations describe. Within one iteration
    /// \\(\Phi(\mathbf{x})\\), \\(A\\), \\(\mathbf{g}\\) and every
    /// trial \\(\Phi(\mathbf{x}+\mathbf{h})\\) — including the ones
    /// spent escalating \\(\mu\\) — must come from ONE objective, or
    /// the gain ratio measures the difference between two different
    /// functions and \\(A\\) is not the Gauss-Newton Hessian of the
    /// cost being compared. Switching models *between* iterations is
    /// ordinary inexact-model refresh and is fully supported here.
    ///
    /// Return `true` to declare that the model just changed. The
    /// driver then re-assembles at `x` under the new model, commits
    /// the result via [`on_step_accepted`](Self::on_step_accepted),
    /// and only then enters the trial loop — so the system it damps
    /// and the trials it compares against are the same objective
    /// again. Note that the objective's *value* may jump across a
    /// refresh: cost monotonicity holds per model, not across one.
    ///
    /// **Cost:** each `true` spends one extra FULL evaluation
    /// (residuals + Jacobian) beyond the trial-driven ones, at the
    /// current point. It is *outside* the
    /// [`max_iterations`](LMConfig::max_iterations) iteration budget —
    /// a refresh does not consume an outer iteration — and is reported
    /// as [`LMSolution::n_model_refreshes`] so a caller sizing a
    /// wall-clock budget can account for it. For a problem whose
    /// assembly is the expensive part (an n-body propagation over an
    /// observation arc), that is the dominant per-refresh cost.
    ///
    /// Return `false` (the default) to keep the current system; no
    /// evaluation is spent.
    fn refresh_model(&mut self, _x: &[f64; N]) -> bool {
        false
    }

    /// Whether the system the driver most recently assembled describes
    /// the problem's **own** objective rather than a stand-in for it.
    ///
    /// Default `true`: an ordinary problem evaluates one function and
    /// every assembly is that function. A problem that answers
    /// [`ResidualProblem::evaluate`] from a cached linearization —
    /// serving \\(f(x_0) + J\,\delta x\\) instead of propagating — returns
    /// `false` for as long as that is so.
    ///
    /// # What the driver does with it
    ///
    /// **An inexact system may not end the solve.** A convergence
    /// criterion firing on one measures the stand-in's stationarity, and
    /// a damping or inner-trial exhaustion on one measures the stand-in's
    /// exhaustion; neither is a statement about the objective. On either,
    /// the driver spends one re-assembly at the current point — the same
    /// evaluation a [`refresh_model`](Self::refresh_model) declaration
    /// buys, counted the same way — and then:
    ///
    /// * re-tests the convergence battery on what came back, latching
    ///   only if a criterion fires there. So every
    ///   \\(\texttt{converged} = \texttt{true}\\) verdict, and the
    ///   \\(\mathbf{x}\\) it certifies, belongs to the objective and not
    ///   to a linearization of it.
    /// * on an exhaustion, discards the damping escalated against the
    ///   stand-in (μ returns to its value at the last committed
    ///   acceptance) and keeps iterating, so a stand-in the driver ran out
    ///   of road on costs a re-anchoring rather than the solve.
    ///
    /// Termination stays guaranteed: a second exhaustion with no step
    /// accepted between the two re-assemblies is the objective's own, and
    /// ends the solve. And every [`LMSolution`] field read off the final
    /// system — the covariance, the gradient norm, the returned
    /// \\(q\\)-norm — is therefore computed on an exact assembly on every
    /// exit path.
    ///
    /// # Contract
    ///
    /// Read immediately after an assembly, and must describe THAT
    /// assembly. A problem that returns `false` must be able to produce an
    /// exact assembly on the next call — the driver's re-assembly is not
    /// optional, and a problem that served a second linearization there
    /// would leave the driver with no way to terminate honestly.
    fn assembly_is_exact(&self) -> bool {
        true
    }

    /// Directional second derivative of the pre-weighted residuals
    /// along `v` at `x`:
    /// \\(\mathbf{r}''_{vv} = d^2\mathbf{r}_w(\mathbf{x} +
    /// t\,\mathbf{v})/dt^2\big|_{t=0}\\), in the SAME row order as
    /// the most recent full evaluation at `x`. Powers geodesic
    /// acceleration ([`LMConfig::geodesic_acceleration`]); a single
    /// one-parameter second-order jet evaluation along `v` suffices —
    /// the full Hessian is never needed.
    ///
    /// `None` (the default) disables acceleration for this step. Must
    /// not commit persistent state (same contract as
    /// [`evaluate_cost`](Self::evaluate_cost)).
    fn second_directional_derivative(&mut self, _x: &[f64; N], _v: &[f64; N]) -> Option<Vec<f64>> {
        None
    }
}

/// Residual-level problem: the driver assembles
/// \\(A = J_w^\top J_w\\), \\(\mathbf{g} = -J_w^\top \mathbf{r}_w\\)
/// from pre-weighted rows and owns the optional prior.
pub trait ResidualProblem<const N: usize>: CostProblem<N> {
    /// Full evaluation (residuals + Jacobian + cost).
    ///
    /// Called at \\(\mathbf{x}_0\\) and at provisionally accepted
    /// points; committed state moves only after it succeeds (a failing
    /// or cost-inconsistent evaluation rolls the acceptance back). The
    /// final COMMITTED evaluation — the one followed by
    /// [`on_step_accepted`](CostProblem::on_step_accepted) —
    /// corresponds to the returned \\(\mathbf{x}\\); a later
    /// `evaluate` call may return `Ok` at a trial point whose
    /// acceptance was then rolled back (`on_step_rejected` fires for
    /// it).
    fn evaluate(&mut self, x: &[f64; N]) -> Result<NLLSEvaluation<N>, Self::Error>;
}

/// Normal-equations-level evaluation: cost, normal matrix, and
/// right-hand side describing **one** objective.
#[derive(Clone, Debug)]
pub struct SystemEvaluation<const N: usize> {
    /// Objective value \\(\Phi\\). Must include problem-composed
    /// priors and, for reduced systems (Schur), be the **profiled**
    /// cost \\(\min_b F(\mathbf{x}, b)\\).
    pub cost: f64,
    /// \\(A\\): symmetric positive semi-definite; the Gauss-Newton
    /// Hessian (×½) of the same objective as `cost`.
    pub normal: [[f64; N]; N],
    /// \\(\mathbf{g}\\): the negated gradient (×½) of the same
    /// objective, i.e. the right-hand side of
    /// \\(A\,\mathbf{h} = \mathbf{g}\\).
    pub rhs: [f64; N],
}

/// Normal-equations-level problem (e.g. a Schur-reduced system).
///
/// Validated on receipt: all entries finite and
/// \\(A_{jj} \ge 0\\), else the solve fails on the
/// [`LMError::InvalidSystem`] axis — a buggy reduction must surface
/// as what it is, not get masked into `DampingExhausted` through a
/// NaN-poisoned scaling diagonal.
pub trait SystemProblem<const N: usize>: CostProblem<N> {
    /// Full evaluation at \\(\mathbf{x}_0\\) and at provisionally
    /// accepted points; same lifecycle contract as
    /// [`ResidualProblem::evaluate`].
    fn evaluate_system(&mut self, x: &[f64; N]) -> Result<SystemEvaluation<N>, Self::Error>;
}

// ── Output ──────────────────────────────────────────────────────────

/// Why the solver stopped.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub enum TerminationReason {
    /// Scaled gradient below `gtol` at an accepted point
    /// (\\(|g_j| \le \texttt{gtol}\, d_j \sqrt{\Phi}\\) for all j).
    GradientTolerance,
    /// The UNDAMPED Gauss-Newton step at the accepted point fell
    /// below `qtol`/`xtol`; \\(\mathbf{x}\\) unchanged. Measured on
    /// the undamped step deliberately: a μ-shrunken step says nothing
    /// about stationarity.
    StepTolerance,
    /// Actual and predicted reduction of the last accepted step both
    /// below `ftol`·Φ with a consistent (ρ ≤ 2) model.
    CostTolerance,
    /// Outer iteration budget exhausted (`converged = false`).
    MaxIterations,
    /// \\(\mu\\) exceeded `mu_max` while rejecting
    /// (`converged = false`).
    DampingExhausted {
        /// The damping value that exceeded the cap.
        mu: f64,
    },
    /// Inner trial budget exhausted without an acceptance
    /// (`converged = false`).
    InnerTrialsExhausted {
        /// Trials consumed in the exhausted iteration.
        trials: usize,
    },
}

/// Why the covariance could not be produced. The solution (`x`, cost,
/// diagnostics) is still returned — a rank-deficient but converged fit
/// keeps its state; it just has no finite σ in the unobservable
/// subspace.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub enum CovarianceFailure {
    /// The (prior-augmented) normal matrix at the returned point is
    /// singular to working precision.
    SingularNormalMatrix,
}

/// Solver output.
///
/// `#[non_exhaustive]`, so struct-literal construction is reserved to
/// the solver — an `LMSolution` built anywhere else would be a
/// fabricated solver output — and future diagnostic fields (like
/// [`n_model_refreshes`](Self::n_model_refreshes), added in 1.11) stay
/// additive for downstream readers instead of breaking exhaustive
/// destructuring silently.
#[derive(Clone, Debug)]
#[non_exhaustive]
pub struct LMSolution<const N: usize> {
    /// The retained iterate (best visited: cost is monotone
    /// non-increasing across accepted steps *of one model* — a
    /// [`refresh_model`](CostProblem::refresh_model) declaration
    /// re-measures the objective, so `cost` may rise across a refresh).
    pub x: [f64; N],
    /// \\((A_{\text{final}})^{-1}\\) at the returned `x`
    /// (data + driver prior). Inversion failure is carried explicitly —
    /// never fabricated as zeros, never discards the solution.
    pub covariance: Result<[[f64; N]; N], CovarianceFailure>,
    /// Full objective \\(\Phi\\) at the returned `x`, from the final
    /// full evaluation (data + prior penalty).
    pub cost: f64,
    /// `cost` minus the driver prior penalty (equals `cost` on the
    /// system path, where the problem owns its priors).
    pub data_cost: f64,
    /// Quadratic form \\(\mathbf{h}^\top A\,\mathbf{h}\\) of the last
    /// **accepted** step under the undamped, prior-augmented \\(A\\)
    /// that produced it — which, after a
    /// [`refresh_model`](CostProblem::refresh_model), is not the
    /// \\(A\\) of the returned system. `None` only when no step was
    /// ever accepted; a refresh does NOT reset it, so `None` always
    /// means what it says.
    ///
    /// A DIAGNOSTIC, not a convergence metric: the accepted step is
    /// μ-damped and possibly clamped, so it is not the quantity any
    /// termination test is decided on and is not comparable to
    /// `qtol`/`xtol`/`ftol` in either direction. Quote
    /// [`final_gn_qnorm`](Self::final_gn_qnorm) when reporting against
    /// a tolerance.
    pub accepted_step_qnorm: Option<f64>,
    /// Quadratic form \\(\mathbf{h}_{\text{GN}}^\top
    /// A\,\mathbf{h}_{\text{GN}}\\) of the **undamped** Gauss-Newton
    /// step at the returned `x` — the quantity the
    /// [`qtol`](LMConfig::qtol) step-convergence test is decided on,
    /// and therefore the only step norm comparable to that tolerance.
    /// `None` when the undamped system is singular at the returned
    /// point, where the test is skipped as well.
    pub final_gn_qnorm: Option<f64>,
    /// \\(\max_j |g_j| / d_j\\) at the returned `x` (diagnostic; the
    /// gradient *test* additionally normalizes by \\(\sqrt{\Phi}\\)).
    pub gradient_norm_scaled: f64,
    /// Outer iterations performed (each holds one accepted system
    /// snapshot; see [`LMConfig::max_iterations`] for the
    /// full-evaluation accounting).
    pub iterations: usize,
    /// Trial-cost evaluations performed.
    pub n_cost_evals: usize,
    /// Trials rejected by the gain-ratio test (finite costs) or by a
    /// failed damped-system factorization (μ raised, no cost
    /// evaluated).
    pub n_rejected_trials: usize,
    /// Trials rejected for `Err`/non-finite costs (including rolled-back
    /// acceptances).
    pub n_invalid_trials: usize,
    /// Trials whose step carried the geodesic-acceleration correction.
    pub n_accelerated_trials: usize,
    /// [`refresh_model`](CostProblem::refresh_model) declarations that
    /// fired, i.e. the number of extra FULL evaluations (residuals +
    /// Jacobian) the driver spent re-assembling at the current point
    /// after the problem swapped its model.
    ///
    /// These are outside the [`max_iterations`](LMConfig::max_iterations)
    /// iteration budget and are counted by no other field, so a caller
    /// sizing a wall-clock budget against assembly cost must add them —
    /// see the evaluation accounting on
    /// [`max_iterations`](LMConfig::max_iterations). Bounded by
    /// [`iterations`](Self::iterations) (a refresh is offered at most
    /// once per outer iteration) and `0` for every problem that leaves
    /// the default hook in place.
    pub n_model_refreshes: usize,
    /// Final damping parameter.
    pub mu_final: f64,
    /// Whether a convergence criterion was met.
    pub converged: bool,
    /// The criterion met, or the budget that ran out.
    pub reason: TerminationReason,
}

// ── Errors ──────────────────────────────────────────────────────────

/// A configuration field violation.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub enum ConfigDefect {
    /// `max_iterations` must be ≥ 1.
    MaxIterationsZero,
    /// `max_inner_trials` must be ≥ 1.
    MaxInnerTrialsZero,
    /// `tau` must be finite and > 0.
    TauNotPositive {
        /// Offending value.
        value: f64,
    },
    /// `mu_max` must be finite and > 0.
    MuMaxNotPositive {
        /// Offending value.
        value: f64,
    },
    /// `tau` must not exceed `mu_max`: \\(\mu_0 \le \tau\\) by
    /// construction, so this guarantees the solve never starts above
    /// the damping cap.
    TauExceedsMuMax {
        /// Configured initial damping scale.
        tau: f64,
        /// Configured damping cap.
        mu_max: f64,
    },
    /// `min_relative_decrease` must be finite and ≥ 0.
    MinRelativeDecreaseNegative {
        /// Offending value.
        value: f64,
    },
    /// `gtol` must be finite and ≥ 0.
    GtolNegative {
        /// Offending value.
        value: f64,
    },
    /// `qtol` must be finite and ≥ 0.
    QtolNegative {
        /// Offending value.
        value: f64,
    },
    /// `xtol` must be finite and ≥ 0.
    XtolNegative {
        /// Offending value.
        value: f64,
    },
    /// `ftol` must be finite and ≥ 0.
    FtolNegative {
        /// Offending value.
        value: f64,
    },
    /// `max_consecutive_invalid` must be ≥ 1.
    MaxConsecutiveInvalidZero,
    /// `avmax` must be finite and > 0.
    AvmaxNotPositive {
        /// Offending value.
        value: f64,
    },
}

/// A prior violation.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub enum PriorDefect {
    /// Non-finite prior mean component.
    NonFiniteMean {
        /// Component index.
        index: usize,
    },
    /// Non-finite inverse-covariance entry.
    NonFiniteCovarianceInv {
        /// Row index.
        row: usize,
        /// Column index.
        col: usize,
    },
    /// Negative inverse-covariance diagonal (not a valid precision
    /// matrix).
    NegativeDiagonal {
        /// Diagonal index.
        index: usize,
        /// Offending value.
        value: f64,
    },
}

/// An invalid value in a full residual-level evaluation.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub enum EvaluationDefect {
    /// `cost` was NaN or infinite.
    Cost,
    /// `cost` was negative — impossible for a sum of squares (plus a
    /// PSD prior penalty). Left unvalidated, a negative cost would
    /// NaN-poison \\(\sqrt{\Phi}\\) in the gradient test and read as
    /// instant false convergence.
    NegativeCost {
        /// Offending value.
        value: f64,
    },
    /// A residual entry was NaN or infinite.
    Residual {
        /// Row index.
        index: usize,
    },
    /// A Jacobian entry was NaN or infinite.
    JacobianEntry {
        /// Row index.
        row: usize,
        /// Column index.
        col: usize,
    },
}

/// An invalid [`SystemEvaluation`].
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub enum SystemDefect {
    /// `cost` was NaN or infinite.
    NonFiniteCost,
    /// A normal-matrix entry was NaN or infinite.
    NonFiniteNormal {
        /// Row index.
        row: usize,
        /// Column index.
        col: usize,
    },
    /// A right-hand-side entry was NaN or infinite.
    NonFiniteRhs {
        /// Component index.
        index: usize,
    },
    /// `cost` was negative — impossible for a real (profiled) least-
    /// squares objective. A Schur consumer whose profiled cost rounds
    /// to a tiny negative through catastrophic cancellation must
    /// handle that deliberately on its side before returning; the
    /// driver never silently clamps.
    NegativeCost {
        /// Offending value.
        value: f64,
    },
    /// A negative normal-matrix diagonal entry — algebraically
    /// impossible for a real Gauss-Newton/Schur reduction; a definite
    /// bug in the problem.
    NegativeDiagonal {
        /// Diagonal index.
        index: usize,
        /// Offending value.
        value: f64,
    },
}

/// Solver error. Variants are split by failure axis; the problem's
/// domain error type `E` is carried through without flattening.
#[derive(Debug)]
#[non_exhaustive]
pub enum LMError<E> {
    /// A configuration field is invalid.
    InvalidConfig {
        /// Which field, and how.
        defect: ConfigDefect,
    },
    /// The supplied prior is invalid.
    InvalidPrior {
        /// Which part, and how.
        defect: PriorDefect,
    },
    /// A [`SystemEvaluation`] failed validation at the initial point
    /// (later structurally-invalid systems with a *negative diagonal*
    /// are also fatal; non-finite systems at provisionally accepted
    /// points roll back instead).
    InvalidSystem {
        /// Outer iteration (0 = initial evaluation).
        iteration: usize,
        /// What was invalid.
        defect: SystemDefect,
    },
    /// A full evaluation returned zero residual rows.
    EmptyResiduals {
        /// Outer iteration (0 = initial evaluation).
        iteration: usize,
    },
    /// Residual and Jacobian row counts disagree.
    DimensionMismatch {
        /// Outer iteration (0 = initial evaluation).
        iteration: usize,
        /// Residual rows returned.
        residuals: usize,
        /// Jacobian rows returned.
        jacobian: usize,
    },
    /// The scaling diagonal is identically zero (zero Jacobian, no
    /// prior): multiplicative damping can never de-singularize it, so
    /// escalating \\(\mu\\) is provably futile.
    ZeroDampingDiagonal {
        /// Outer iteration (0 = initial evaluation).
        iteration: usize,
    },
    /// The problem failed at \\(\mathbf{x}_0\\) — there is no valid
    /// reference point to retreat to.
    InitialEvaluationFailed {
        /// The domain error.
        source: E,
    },
    /// The re-assembly that a
    /// [`refresh_model`](CostProblem::refresh_model) declaration
    /// requires failed at the accepted point. There is no retreat: the
    /// problem has already moved to the new model, so the system the
    /// driver still holds describes an objective the problem no longer
    /// evaluates.
    ModelRefreshFailed {
        /// Outer iteration whose refresh failed.
        iteration: usize,
        /// The domain error.
        source: E,
    },
    /// The re-assembly forced by an inexact system came back inexact
    /// too, so the driver has no exact system to judge and no way to
    /// terminate honestly. See
    /// [`assembly_is_exact`](CostProblem::assembly_is_exact) — a problem
    /// that declares an assembly inexact must be able to produce an exact
    /// one when the driver asks.
    InexactAssemblyPersisted {
        /// Outer iteration whose forced re-assembly stayed inexact.
        iteration: usize,
    },
    /// Non-finite value in a full evaluation the driver cannot retreat
    /// from: the **initial** one, or the re-assembly a
    /// [`refresh_model`](CostProblem::refresh_model) declaration
    /// requires (there the problem has already switched models, so
    /// there is no earlier system left to fall back to). A non-finite
    /// TRIAL evaluation rolls the acceptance back instead of surfacing
    /// here.
    InvalidEvaluation {
        /// Outer iteration (0 = initial evaluation).
        iteration: usize,
        /// What was non-finite.
        defect: EvaluationDefect,
    },
    /// [`LMConfig::square_root_solve`] was requested, the problem is on
    /// the residual path, and the assembled system carried no
    /// least-squares rows.
    ///
    /// The only way a residual-path assembly withholds its rows is a
    /// prior whose precision matrix cannot be written as \\(S^\top S\\)
    /// — one with a genuinely negative eigenvalue, which is not a
    /// precision matrix at all. Reported rather than absorbed: without
    /// this the
    /// step solve returns `None` on every trial, the driver reads that
    /// as a rejection, and the solve terminates with
    /// [`TerminationReason::DampingExhausted`] at enormous \\(\mu\\)
    /// having never once evaluated the objective — a diagnosis that
    /// points at the damping policy instead of at the prior.
    SquareRootSystemUnavailable {
        /// Outer iteration that could not build a step.
        iteration: usize,
    },
    /// `max_consecutive_invalid` consecutive trial evaluations failed
    /// (`Err`, non-finite, or rolled-back acceptances). The best
    /// visited (= current accepted) state is carried for triage.
    PersistentInvalidTrials {
        /// Outer iteration of the last failure.
        iteration: usize,
        /// Consecutive failures observed.
        consecutive: usize,
        /// The last domain error, if the failures carried one.
        last_source: Option<E>,
        /// Best (= last accepted) iterate.
        best_x: Vec<f64>,
        /// Objective at `best_x`.
        best_cost: f64,
    },
}

impl<E: std::fmt::Display> std::fmt::Display for LMError<E> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::InvalidConfig { defect } => write!(f, "invalid configuration: {defect:?}"),
            Self::InvalidPrior { defect } => write!(f, "invalid prior: {defect:?}"),
            Self::InvalidSystem { iteration, defect } => {
                write!(
                    f,
                    "invalid system evaluation at iteration {iteration}: {defect:?}"
                )
            }
            Self::EmptyResiduals { iteration } => {
                write!(
                    f,
                    "evaluation returned no residuals at iteration {iteration}"
                )
            }
            Self::DimensionMismatch {
                iteration,
                residuals,
                jacobian,
            } => write!(
                f,
                "dimension mismatch at iteration {iteration}: {residuals} residuals vs \
                 {jacobian} Jacobian rows"
            ),
            Self::ZeroDampingDiagonal { iteration } => write!(
                f,
                "scaling diagonal identically zero at iteration {iteration} (zero Jacobian, \
                 no prior): damping cannot regularize this system"
            ),
            Self::InitialEvaluationFailed { source } => {
                write!(f, "initial evaluation failed: {source}")
            }
            Self::ModelRefreshFailed { iteration, source } => write!(
                f,
                "re-assembly after a model refresh failed at iteration {iteration}: {source}"
            ),
            Self::InexactAssemblyPersisted { iteration } => write!(
                f,
                "the re-assembly forced at iteration {iteration} came back inexact: a problem \
                 that declares an assembly inexact must be able to produce an exact one when \
                 asked, or no verdict the driver reaches describes its objective"
            ),
            Self::InvalidEvaluation { iteration, defect } => {
                write!(f, "invalid evaluation at iteration {iteration}: {defect:?}")
            }
            Self::SquareRootSystemUnavailable { iteration } => write!(
                f,
                "square-root solve requested but the assembled system carried no least-squares \
                 rows at iteration {iteration}: the prior's precision matrix has no square root \
                 S with S^T S = P0^-1, so the prior cannot enter the factorization"
            ),
            Self::PersistentInvalidTrials {
                iteration,
                consecutive,
                best_cost,
                ..
            } => write!(
                f,
                "{consecutive} consecutive invalid trial evaluations at iteration {iteration} \
                 (best retained cost {best_cost})"
            ),
        }
    }
}

impl<E: std::error::Error + 'static> std::error::Error for LMError<E> {
    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
        match self {
            Self::InitialEvaluationFailed { source } => Some(source),
            Self::ModelRefreshFailed { source, .. } => Some(source),
            Self::PersistentInvalidTrials {
                last_source: Some(source),
                ..
            } => Some(source),
            _ => None,
        }
    }
}

// ── Internal: assembled system + sources ────────────────────────────

/// One full evaluation, reduced to the common driver currency.
struct AssembledSystem<const N: usize> {
    /// Full objective Φ (prior penalty included on the residual path).
    cost: f64,
    /// Objective minus the driver prior penalty.
    data_cost: f64,
    /// Symmetrized, prior-augmented normal matrix A.
    normal: [[f64; N]; N],
    /// Prior-augmented right-hand side g (so that A h = g).
    rhs: [f64; N],
    /// Pre-weighted Jacobian rows (residual path only) — needed for
    /// the geodesic-acceleration product \\(J^\top \mathbf{r}''_{vv}\\).
    /// `None` on the normal-equations path, which therefore cannot
    /// use acceleration.
    jacobian: Option<Vec<[f64; N]>>,
    /// The least-squares system in SQUARE-ROOT form: stacked rows
    /// \([J; L^\top]\) with right-hand side \([-r; -L^\top(x-m)]\),
    /// where \(P_0^{-1} = L L^\top\) contributes the prior.
    ///
    /// Carried alongside `normal`/`rhs` rather than replacing them: the
    /// gain ratio, the gradient test and the predicted reduction are all
    /// defined on \(A\) and \(g\), so the square-root path changes only
    /// how the STEP is computed. `None` on the normal-equations path,
    /// which never sees rows.
    ls_rows: Option<(Vec<[f64; N]>, Vec<f64>)>,
}

/// How a full evaluation failed.
enum AssembleFailure<E> {
    /// Problem-domain error: fatal at x0, rolls back elsewhere.
    Domain(E),
    /// Non-finite residual-path value: fatal at x0, rolls back
    /// elsewhere.
    NonFinite(EvaluationDefect),
    /// Non-finite system-path value: fatal at x0, rolls back elsewhere.
    SystemNonFinite(SystemDefect),
    /// Contract violation (empty/mismatched rows, negative diagonal):
    /// always fatal.
    Hard(HardDefect),
}

enum HardDefect {
    EmptyResiduals,
    DimensionMismatch { residuals: usize, jacobian: usize },
    NegativeDiagonal { index: usize, value: f64 },
}

impl HardDefect {
    fn into_error<E>(self, iteration: usize) -> LMError<E> {
        match self {
            Self::EmptyResiduals => LMError::EmptyResiduals { iteration },
            Self::DimensionMismatch {
                residuals,
                jacobian,
            } => LMError::DimensionMismatch {
                iteration,
                residuals,
                jacobian,
            },
            Self::NegativeDiagonal { index, value } => LMError::InvalidSystem {
                iteration,
                defect: SystemDefect::NegativeDiagonal { index, value },
            },
        }
    }
}

/// Adapter unifying the residual and system paths for the driver core.
trait SystemSource<const N: usize> {
    type Error: std::error::Error + 'static;

    fn assemble(
        &mut self,
        x: &[f64; N],
    ) -> Result<AssembledSystem<N>, AssembleFailure<Self::Error>>;
    fn trial_cost(&mut self, x: &[f64; N]) -> Result<f64, Self::Error>;
    fn constrain(&mut self, x: &[f64; N], delta: &mut [f64; N]);
    fn accepted(&mut self, x: &[f64; N]);
    fn rejected(&mut self, x_trial: &[f64; N]);
    fn refresh_model(&mut self, x: &[f64; N]) -> bool;
    fn assembly_is_exact(&self) -> bool;
    fn second_directional_derivative(&mut self, x: &[f64; N], v: &[f64; N]) -> Option<Vec<f64>>;
}

/// Residual-level adapter: assembles normal equations, owns the prior
/// (an exactly-symmetrized copy — see [`solve`]).
struct ResidualSource<'a, P, const N: usize> {
    problem: &'a mut P,
    prior: Option<NLLSPrior<N>>,
}

impl<'a, P: ResidualProblem<N>, const N: usize> ResidualSource<'a, P, N> {
    /// Prior penalty \\((\mathbf{x}-\mathbf{m})^\top P_0^{-1} (\mathbf{x}-\mathbf{m})\\),
    /// fixed evaluation order.
    fn prior_penalty(&self, x: &[f64; N]) -> f64 {
        let Some(p) = &self.prior else { return 0.0 };
        let mut penalty = 0.0_f64;
        for i in 0..N {
            let mut row = 0.0_f64;
            for j in 0..N {
                row += p.covariance_inv[i][j] * (x[j] - p.mean[j]);
            }
            penalty += (x[i] - p.mean[i]) * row;
        }
        penalty
    }
}

impl<'a, P: ResidualProblem<N>, const N: usize> SystemSource<N> for ResidualSource<'a, P, N> {
    type Error = P::Error;

    fn assemble(
        &mut self,
        x: &[f64; N],
    ) -> Result<AssembledSystem<N>, AssembleFailure<Self::Error>> {
        let eval = self.problem.evaluate(x).map_err(AssembleFailure::Domain)?;

        if eval.residuals.is_empty() {
            return Err(AssembleFailure::Hard(HardDefect::EmptyResiduals));
        }
        if eval.residuals.len() != eval.jacobian.len() {
            return Err(AssembleFailure::Hard(HardDefect::DimensionMismatch {
                residuals: eval.residuals.len(),
                jacobian: eval.jacobian.len(),
            }));
        }
        if !eval.cost.is_finite() {
            return Err(AssembleFailure::NonFinite(EvaluationDefect::Cost));
        }
        if eval.cost < 0.0 {
            return Err(AssembleFailure::NonFinite(EvaluationDefect::NegativeCost {
                value: eval.cost,
            }));
        }
        for (i, r) in eval.residuals.iter().enumerate() {
            if !r.is_finite() {
                return Err(AssembleFailure::NonFinite(EvaluationDefect::Residual {
                    index: i,
                }));
            }
        }
        for (i, row) in eval.jacobian.iter().enumerate() {
            for (j, v) in row.iter().enumerate() {
                if !v.is_finite() {
                    return Err(AssembleFailure::NonFinite(
                        EvaluationDefect::JacobianEntry { row: i, col: j },
                    ));
                }
            }
        }

        // Normal equations: A = JᵀJ, g = -Jᵀr (fixed accumulation order).
        let mut normal = [[0.0_f64; N]; N];
        let mut rhs = [0.0_f64; N];
        for (r_i, j_i) in eval.residuals.iter().zip(eval.jacobian.iter()) {
            for j in 0..N {
                for k in 0..N {
                    normal[j][k] += j_i[j] * j_i[k];
                }
                rhs[j] -= j_i[j] * r_i;
            }
        }
        let mut normal = mat_symmetrize(&normal);

        // Prior augmentation: A += P0⁻¹, g -= P0⁻¹(x − m). The prior
        // was symmetrized on receipt, so `normal` stays exactly
        // symmetric and the solved system agrees with the quadratic
        // model used for the gain ratio.
        if let Some(p) = &self.prior {
            for i in 0..N {
                for j in 0..N {
                    normal[i][j] += p.covariance_inv[i][j];
                }
                let mut delta_weighted = 0.0_f64;
                for j in 0..N {
                    delta_weighted += p.covariance_inv[i][j] * (x[j] - p.mean[j]);
                }
                rhs[i] -= delta_weighted;
            }
        }

        let data_cost = eval.cost;
        let cost = data_cost + self.prior_penalty(x);
        if !cost.is_finite() {
            return Err(AssembleFailure::NonFinite(EvaluationDefect::Cost));
        }
        if cost < 0.0 {
            return Err(AssembleFailure::NonFinite(EvaluationDefect::NegativeCost {
                value: cost,
            }));
        }

        // Square-root form of the SAME system: data rows as given, plus
        // the prior as appended rows. A Cholesky of P0^-1 = L L^T makes
        // the prior penalty (x-m)^T P0^-1 (x-m) equal to ||L^T (x-m)||^2,
        // so L^T rows with residual L^T(x-m) reproduce it exactly — the
        // prior enters the factorization instead of being added to a
        // squared matrix.
        let ls_rows = {
            let mut rows = eval.jacobian.clone();
            let mut b: Vec<f64> = eval.residuals.iter().map(|r| -r).collect();
            let mut ok = true;
            if let Some(p) = &self.prior {
                match prior_rows(&p.covariance_inv) {
                    Some(factor) => {
                        for row in factor {
                            let resid: f64 = (0..N).map(|j| row[j] * (x[j] - p.mean[j])).sum();
                            rows.push(row);
                            b.push(-resid);
                        }
                    }
                    // A precision matrix with a genuinely negative
                    // eigenvalue is not a sum of squares in ANY basis,
                    // so no set of rows reproduces its penalty.
                    // Withhold rather than silently drop the prior from
                    // the step; `solve_core` surfaces the withholding.
                    None => ok = false,
                }
            }
            ok.then_some((rows, b))
        };

        Ok(AssembledSystem {
            cost,
            data_cost,
            normal,
            rhs,
            jacobian: Some(eval.jacobian),
            ls_rows,
        })
    }

    fn trial_cost(&mut self, x: &[f64; N]) -> Result<f64, Self::Error> {
        let data = self.problem.evaluate_cost(x)?;
        // A non-finite data cost propagates as non-finite total: the
        // driver rejects it (never substitutes).
        Ok(data + self.prior_penalty(x))
    }

    fn constrain(&mut self, x: &[f64; N], delta: &mut [f64; N]) {
        self.problem.constrain_step(x, delta);
    }

    fn accepted(&mut self, x: &[f64; N]) {
        self.problem.on_step_accepted(x);
    }

    fn rejected(&mut self, x_trial: &[f64; N]) {
        self.problem.on_step_rejected(x_trial);
    }

    fn refresh_model(&mut self, x: &[f64; N]) -> bool {
        self.problem.refresh_model(x)
    }

    fn assembly_is_exact(&self) -> bool {
        self.problem.assembly_is_exact()
    }

    fn second_directional_derivative(&mut self, x: &[f64; N], v: &[f64; N]) -> Option<Vec<f64>> {
        self.problem.second_directional_derivative(x, v)
    }
}

/// Normal-equations-level adapter: the problem owns its objective
/// (priors included); the driver only validates.
struct DirectSource<'a, P, const N: usize> {
    problem: &'a mut P,
}

impl<'a, P: SystemProblem<N>, const N: usize> SystemSource<N> for DirectSource<'a, P, N> {
    type Error = P::Error;

    fn assemble(
        &mut self,
        x: &[f64; N],
    ) -> Result<AssembledSystem<N>, AssembleFailure<Self::Error>> {
        let sys = self
            .problem
            .evaluate_system(x)
            .map_err(AssembleFailure::Domain)?;

        if !sys.cost.is_finite() {
            return Err(AssembleFailure::SystemNonFinite(
                SystemDefect::NonFiniteCost,
            ));
        }
        if sys.cost < 0.0 {
            return Err(AssembleFailure::SystemNonFinite(
                SystemDefect::NegativeCost { value: sys.cost },
            ));
        }
        for i in 0..N {
            for j in 0..N {
                if !sys.normal[i][j].is_finite() {
                    return Err(AssembleFailure::SystemNonFinite(
                        SystemDefect::NonFiniteNormal { row: i, col: j },
                    ));
                }
            }
            if !sys.rhs[i].is_finite() {
                return Err(AssembleFailure::SystemNonFinite(
                    SystemDefect::NonFiniteRhs { index: i },
                ));
            }
        }
        for i in 0..N {
            if sys.normal[i][i] < 0.0 {
                return Err(AssembleFailure::Hard(HardDefect::NegativeDiagonal {
                    index: i,
                    value: sys.normal[i][i],
                }));
            }
        }

        Ok(AssembledSystem {
            cost: sys.cost,
            data_cost: sys.cost,
            normal: mat_symmetrize(&sys.normal),
            rhs: sys.rhs,
            jacobian: None,
            ls_rows: None,
        })
    }

    fn trial_cost(&mut self, x: &[f64; N]) -> Result<f64, Self::Error> {
        self.problem.evaluate_cost(x)
    }

    fn constrain(&mut self, x: &[f64; N], delta: &mut [f64; N]) {
        self.problem.constrain_step(x, delta);
    }

    fn accepted(&mut self, x: &[f64; N]) {
        self.problem.on_step_accepted(x);
    }

    fn rejected(&mut self, x_trial: &[f64; N]) {
        self.problem.on_step_rejected(x_trial);
    }

    fn refresh_model(&mut self, x: &[f64; N]) -> bool {
        self.problem.refresh_model(x)
    }

    fn assembly_is_exact(&self) -> bool {
        self.problem.assembly_is_exact()
    }

    fn second_directional_derivative(&mut self, x: &[f64; N], v: &[f64; N]) -> Option<Vec<f64>> {
        self.problem.second_directional_derivative(x, v)
    }
}

// ── Public entry points ─────────────────────────────────────────────

/// Solve a residual-level nonlinear least-squares problem.
///
/// The optional `prior` is driver-owned: it augments the normal
/// equations **and** enters both compared costs, so the acceptance
/// test judges the same MAP objective the step minimizes.
pub fn solve<P: ResidualProblem<N>, const N: usize>(
    problem: &mut P,
    x0: [f64; N],
    config: &LMConfig,
    prior: Option<&NLLSPrior<N>>,
) -> Result<LMSolution<N>, LMError<P::Error>> {
    validate_config(config)?;
    let prior = match prior {
        Some(p) => {
            validate_prior(p)?;
            // Symmetrize the precision matrix once on receipt
            // (bitwise no-op for exactly symmetric input): the solved
            // system, the gain-ratio model, and the penalty must all
            // see the SAME quadratic form, and real priors derived
            // from an LU-based inversion are asymmetric at the ulp
            // level.
            Some(NLLSPrior {
                mean: p.mean,
                covariance_inv: mat_symmetrize(&p.covariance_inv),
            })
        }
        None => None,
    };
    let mut source = ResidualSource { problem, prior };
    solve_core(&mut source, x0, config)
}

/// Solve a normal-equations-level problem (e.g. Schur-reduced).
///
/// No prior argument: the problem owns its complete objective,
/// priors composed in.
pub fn solve_system<P: SystemProblem<N>, const N: usize>(
    problem: &mut P,
    x0: [f64; N],
    config: &LMConfig,
) -> Result<LMSolution<N>, LMError<P::Error>> {
    validate_config(config)?;
    let mut source = DirectSource { problem };
    solve_core(&mut source, x0, config)
}

/// Closure convenience for [`solve`]. Both closures are required —
/// a defaulted trial cost would be a hidden full-evaluation fallback.
pub fn solve_nlls<E, FEval, FCost, const N: usize>(
    eval: FEval,
    cost: FCost,
    x0: [f64; N],
    config: &LMConfig,
    prior: Option<&NLLSPrior<N>>,
) -> Result<LMSolution<N>, LMError<E>>
where
    E: std::error::Error + 'static,
    FEval: FnMut(&[f64; N]) -> Result<NLLSEvaluation<N>, E>,
    FCost: FnMut(&[f64; N]) -> Result<f64, E>,
{
    struct ClosureProblem<FEval, FCost> {
        eval: FEval,
        cost: FCost,
    }
    impl<E, FEval, FCost, const N: usize> CostProblem<N> for ClosureProblem<FEval, FCost>
    where
        E: std::error::Error + 'static,
        FEval: FnMut(&[f64; N]) -> Result<NLLSEvaluation<N>, E>,
        FCost: FnMut(&[f64; N]) -> Result<f64, E>,
    {
        type Error = E;
        fn evaluate_cost(&mut self, x: &[f64; N]) -> Result<f64, E> {
            (self.cost)(x)
        }
    }
    impl<E, FEval, FCost, const N: usize> ResidualProblem<N> for ClosureProblem<FEval, FCost>
    where
        E: std::error::Error + 'static,
        FEval: FnMut(&[f64; N]) -> Result<NLLSEvaluation<N>, E>,
        FCost: FnMut(&[f64; N]) -> Result<f64, E>,
    {
        fn evaluate(&mut self, x: &[f64; N]) -> Result<NLLSEvaluation<N>, E> {
            (self.eval)(x)
        }
    }
    let mut problem = ClosureProblem { eval, cost };
    solve(&mut problem, x0, config, prior)
}

// ── Validation ──────────────────────────────────────────────────────

fn validate_config<E>(config: &LMConfig) -> Result<(), LMError<E>> {
    let defect = if config.max_iterations == 0 {
        Some(ConfigDefect::MaxIterationsZero)
    } else if config.max_inner_trials == 0 {
        Some(ConfigDefect::MaxInnerTrialsZero)
    } else if !(config.tau.is_finite() && config.tau > 0.0) {
        Some(ConfigDefect::TauNotPositive { value: config.tau })
    } else if !(config.mu_max.is_finite() && config.mu_max > 0.0) {
        Some(ConfigDefect::MuMaxNotPositive {
            value: config.mu_max,
        })
    } else if config.tau > config.mu_max {
        // μ₀ ≤ τ by construction (iteration-0 scaling), so τ ≤ mu_max
        // guarantees the solve never STARTS above the damping cap.
        Some(ConfigDefect::TauExceedsMuMax {
            tau: config.tau,
            mu_max: config.mu_max,
        })
    } else if !(config.min_relative_decrease.is_finite() && config.min_relative_decrease >= 0.0) {
        Some(ConfigDefect::MinRelativeDecreaseNegative {
            value: config.min_relative_decrease,
        })
    } else if !(config.gtol.is_finite() && config.gtol >= 0.0) {
        Some(ConfigDefect::GtolNegative { value: config.gtol })
    } else if !(config.qtol.is_finite() && config.qtol >= 0.0) {
        Some(ConfigDefect::QtolNegative { value: config.qtol })
    } else if !(config.xtol.is_finite() && config.xtol >= 0.0) {
        Some(ConfigDefect::XtolNegative { value: config.xtol })
    } else if !(config.ftol.is_finite() && config.ftol >= 0.0) {
        Some(ConfigDefect::FtolNegative { value: config.ftol })
    } else if config.max_consecutive_invalid == 0 {
        Some(ConfigDefect::MaxConsecutiveInvalidZero)
    } else if !(config.avmax.is_finite() && config.avmax > 0.0) {
        Some(ConfigDefect::AvmaxNotPositive {
            value: config.avmax,
        })
    } else {
        None
    };
    match defect {
        Some(defect) => Err(LMError::InvalidConfig { defect }),
        None => Ok(()),
    }
}

fn validate_prior<E, const N: usize>(prior: &NLLSPrior<N>) -> Result<(), LMError<E>> {
    for (i, m) in prior.mean.iter().enumerate() {
        if !m.is_finite() {
            return Err(LMError::InvalidPrior {
                defect: PriorDefect::NonFiniteMean { index: i },
            });
        }
    }
    for i in 0..N {
        for j in 0..N {
            if !prior.covariance_inv[i][j].is_finite() {
                return Err(LMError::InvalidPrior {
                    defect: PriorDefect::NonFiniteCovarianceInv { row: i, col: j },
                });
            }
        }
        if prior.covariance_inv[i][i] < 0.0 {
            return Err(LMError::InvalidPrior {
                defect: PriorDefect::NegativeDiagonal {
                    index: i,
                    value: prior.covariance_inv[i][i],
                },
            });
        }
    }
    Ok(())
}

/// Rows \\(S\\) with \\(S^\top S = P_0^{-1}\\): the prior's contribution
/// to the square-root system, so the penalty
/// \\((\mathbf{x}-\mathbf{m})^\top P_0^{-1}(\mathbf{x}-\mathbf{m})\\)
/// is reproduced exactly as \\(\lVert S(\mathbf{x}-\mathbf{m})\rVert^2\\)
/// without the precision matrix ever being added to a squared matrix.
///
/// # Why this is not just a Cholesky
///
/// A **partial** prior — information on some components and none on the
/// others — is a perfectly ordinary thing to want (a prior on position
/// but not velocity, or on the non-gravitational parameters alone), and
/// [`validate_prior`] admits it: a zero precision row means "no
/// information", not "invalid". Its precision matrix is positive
/// SEMI-definite, and a semi-definite matrix has no Cholesky factor —
/// the diagonal hits zero and the factorization divides by it.
///
/// The normal-equations path never noticed, because it only ever ADDS
/// \\(P_0^{-1}\\) to \\(A\\). The square-root path needs a factor, and a
/// Cholesky-only construction silently withheld the entire square-root
/// system for such a prior: every step solve then returned `None`, the
/// driver read that as a trial rejection, and μ escalated to
/// `DampingExhausted` without the objective ever being evaluated once
/// (empyrean-wv13f).
///
/// The symmetric eigendecomposition has no such gap:
/// \\(P_0^{-1} = V\Lambda V^\top\\) gives rows
/// \\(\sqrt{\lambda_k}\,\mathbf{v}_k^\top\\) over the positive
/// eigenvalues, which reproduce the penalty exactly and simply contribute
/// nothing along the unconstrained directions. Cholesky is still tried
/// first: it is cheaper, and it keeps every positive-definite prior
/// bit-for-bit on the rows it already produced.
///
/// Every positive eigenvalue contributes a row, however small: a weakly
/// informative direction is still information, and dropping it below some
/// threshold would quietly discard part of the prior. Only an exact zero
/// contributes nothing. `None` is reserved for a precision matrix with a
/// genuinely negative eigenvalue — one more than [`PRIOR_EIG_REL_TOL`]
/// below zero relative to the largest, so the rounding dust of an
/// inverted covariance does not trip it. Such a matrix is not a precision
/// matrix at all and is not a sum of squares in any basis, so no set of
/// rows can represent it.
fn prior_rows<const N: usize>(covariance_inv: &[[f64; N]; N]) -> Option<Vec<[f64; N]>> {
    if let Some(l) = mat_cholesky(covariance_inv) {
        // Row i of L^T is column i of L.
        return Some(
            (0..N)
                .map(|i| {
                    let mut row = [0.0_f64; N];
                    for (j, slot) in row.iter_mut().enumerate() {
                        *slot = l[j][i];
                    }
                    row
                })
                .collect(),
        );
    }

    let (eigenvalues, vectors) = mat_symmetric_eigen(covariance_inv)?;
    let scale = eigenvalues.iter().fold(0.0_f64, |m, e| m.max(e.abs()));
    let floor = -PRIOR_EIG_REL_TOL * scale;
    let mut rows = Vec::with_capacity(N);
    for (k, &lambda) in eigenvalues.iter().enumerate() {
        if lambda < floor {
            return None;
        }
        if lambda <= 0.0 {
            continue;
        }
        let root = lambda.sqrt();
        let mut row = [0.0_f64; N];
        for (j, slot) in row.iter_mut().enumerate() {
            *slot = root * vectors[j][k];
        }
        rows.push(row);
    }
    Some(rows)
}

/// Relative floor below which a negative eigenvalue of a precision
/// matrix is rounding dust from inverting a covariance rather than a
/// genuinely indefinite prior.
const PRIOR_EIG_REL_TOL: f64 = 1e-12;

// ── Numerics helpers (fixed-order, deterministic) ───────────────────

/// Update the Moré scaling diagonal in place:
/// \\(d_j \leftarrow \max(d_j, \sqrt{A_{jj}})\\) (running max, never
/// decreasing). The stored scales stay HONEST — no floor is baked in,
/// so genuinely small scales (e.g. ~1e-14 non-gravitational
/// parameters next to O(1) positions) keep their true magnitudes in
/// every norm and test. The floor is applied only at division sites
/// via [`effective_scale`]. Returns the maximum diagonal, which is 0
/// iff the system is identically unscalable.
fn update_scaling<const N: usize>(d: &mut [f64; N], normal: &[[f64; N]; N]) -> f64 {
    let mut d_max = 0.0_f64;
    for j in 0..N {
        let col = normal[j][j].sqrt();
        if col > d[j] {
            d[j] = col;
        }
        if d[j] > d_max {
            d_max = d[j];
        }
    }
    d_max
}

/// Scale used where a division by \\(d_j\\) must be protected from an
/// exactly-zero column: honest \\(d_j\\) when positive, else the
/// relative floor \\(10^{-12}\max_k d_k\\). The floor exists ONLY to
/// keep μD² damping and the equilibration well-defined on zero
/// columns; it never overrides an honest nonzero scale.
#[inline]
fn effective_scale(d_j: f64, d_max: f64) -> f64 {
    if d_j > 0.0 { d_j } else { D_FLOOR_REL * d_max }
}

/// The damped solves run in the **equilibrated basis**
/// \\((B + \mu I)\,\mathbf{h}' = D^{-1}\mathbf{g}\\),
/// \\(B = D^{-1} A D^{-1}\\) — algebraically identical to
/// \\((A + \mu D^2)\,\mathbf{h} = \mathbf{g}\\), but with an O(1)
/// diagonal even when the parameters mix units spanning many orders
/// of magnitude.
///
/// Cholesky factor of the equilibrated damped system
/// \\(B + \mu I = D^{-1} A D^{-1} + \mu I\\). One factorization
/// serves both the velocity and (geodesic) acceleration solves.
fn damped_factor<const N: usize>(
    normal: &[[f64; N]; N],
    d: &[f64; N],
    d_max: f64,
    mu: f64,
) -> Option<[[f64; N]; N]> {
    let mut b = [[0.0_f64; N]; N];
    for i in 0..N {
        let di = effective_scale(d[i], d_max);
        for j in 0..N {
            b[i][j] = normal[i][j] / (di * effective_scale(d[j], d_max));
            // mat_cholesky has no NaN guard; a subnormal-scale 0/0
            // corner must surface here as an ordinary solve failure.
            if !b[i][j].is_finite() {
                return None;
            }
        }
        b[i][i] += mu;
    }
    mat_cholesky(&b)
}

/// Solve with a previously computed [`damped_factor`]:
/// \\((A + \mu D^2)\,\mathbf{h} = \mathbf{g}\\) through the
/// equilibrated triangular solves, mapping back via
/// \\(\mathbf{h} = D^{-1}\mathbf{h}'\\).
fn solve_with_factor<const N: usize>(
    l: &[[f64; N]; N],
    rhs: &[f64; N],
    d: &[f64; N],
    d_max: f64,
) -> Option<[f64; N]> {
    let mut scaled_rhs = [0.0_f64; N];
    for i in 0..N {
        scaled_rhs[i] = rhs[i] / effective_scale(d[i], d_max);
    }

    // Forward solve L y = D⁻¹g.
    let mut y = [0.0_f64; N];
    for i in 0..N {
        let mut sum = scaled_rhs[i];
        for k in 0..i {
            sum -= l[i][k] * y[k];
        }
        y[i] = sum / l[i][i];
    }
    // Back solve Lᵀ h' = y.
    let mut hp = [0.0_f64; N];
    for i in (0..N).rev() {
        let mut sum = y[i];
        for k in (i + 1)..N {
            sum -= l[k][i] * hp[k];
        }
        hp[i] = sum / l[i][i];
    }

    let mut h = [0.0_f64; N];
    for i in 0..N {
        h[i] = hp[i] / effective_scale(d[i], d_max);
        if !h[i].is_finite() {
            return None;
        }
    }
    Some(h)
}

/// Solve \\((A + \mu D^2)\,\mathbf{h} = \mathbf{g}\\) —
/// factorization plus triangular solves. Returns `None` when the
/// factorization fails or the solution is non-finite; the driver
/// treats that as a trial rejection (raise \\(\mu\\), retry) —
/// never as a fallback to another solver.
fn solve_damped<const N: usize>(
    normal: &[[f64; N]; N],
    rhs: &[f64; N],
    d: &[f64; N],
    d_max: f64,
    mu: f64,
) -> Option<[f64; N]> {
    let l = damped_factor(normal, d, d_max, mu)?;
    solve_with_factor(&l, rhs, d, d_max)
}

/// Predicted reduction of the objective for the **actual** (possibly
/// clamped) step, from the general quadratic model:
/// \\(\mathrm{pred} = 2\,\mathbf{h}^\top\mathbf{g} - \mathbf{h}^\top A\,\mathbf{h}\\)
/// (the factor 2 reflects the unhalved-cost convention
/// \\(\Phi = \sum r^2\\)). Valid for arbitrary steps; the exact-step
/// shortcut \\(\mathbf{h}^\top(\mu D^2 \mathbf{h} + \mathbf{g})\\) is
/// not used because clamps invalidate it.
/// Predicted reduction computed from the least-squares ROWS rather than from
/// \\(A\\) and \\(g\\).
///
/// Algebraically identical to [`predicted_reduction`]:
/// \\(2h^\top g - h^\top A h = -2 r^\top J h - \lVert Jh \rVert^2\\)
/// when \\(A = J^\top J\\) and \\(g = -J^\top r\\). Numerically it is not the
/// same thing at all — it never touches the squared matrix, so it stays
/// meaningful at a condition number where \\(h^\top A h\\) has lost every
/// significant digit.
///
/// This matters because the gain ratio decides acceptance. A step computed
/// stably against a predicted reduction computed unstably is rejected on
/// noise, and the damping escalates until it is exhausted — which is the
/// failure this exists to remove.
///
/// `rows` and `b` are the stacked system with \\(b = -r\\).
/// The exact second-order correction to a predicted reduction:
/// \\(h^\top S h = \sum_i r_i (h^\top \nabla^2 r_i h)\\), with
/// \\(S = \sum_i r_i \nabla^2 r_i\\) the part of the least-squares Hessian
/// Gauss-Newton drops.
///
/// `w` is the second directional derivative of the residual vector along
/// \\(h\\) — \\((h^\top \nabla^2 r_i h)_i\\) — and `b` is the least-squares
/// right-hand side, \\(b = -r\\). So
/// \\(h^\top S h = \sum_i r_i w_i = -\sum_i b_i w_i\\), and the corrected
/// model reduction is \\(\text{pred} - h^\top S h = \text{pred} + \sum_i b_i w_i\\).
///
/// `None` on a length mismatch — a partial correction would be worse than
/// none, since it would silently apply curvature from some residuals and
/// not others.
fn second_order_correction(w: &[f64], b: &[f64]) -> Option<f64> {
    if w.len() != b.len() {
        return None;
    }
    Some(w.iter().zip(b).map(|(w_i, b_i)| w_i * b_i).sum())
}

fn predicted_reduction_rows<const N: usize>(h: &[f64; N], rows: &[[f64; N]], b: &[f64]) -> f64 {
    let mut acc = 0.0_f64;
    for (row, b_i) in rows.iter().zip(b) {
        let mut jh = 0.0_f64;
        for j in 0..N {
            jh += row[j] * h[j];
        }
        // b = -r, so -2 r (Jh) = 2 b (Jh).
        acc += 2.0 * b_i * jh - jh * jh;
    }
    acc
}

fn predicted_reduction<const N: usize>(
    h: &[f64; N],
    normal: &[[f64; N]; N],
    rhs: &[f64; N],
) -> f64 {
    let mut hg = 0.0_f64;
    let mut hah = 0.0_f64;
    for i in 0..N {
        hg += h[i] * rhs[i];
        let mut row = 0.0_f64;
        for j in 0..N {
            row += normal[i][j] * h[j];
        }
        hah += h[i] * row;
    }
    2.0 * hg - hah
}

/// Quadratic form \\(\mathbf{h}^\top A\,\mathbf{h}\\) (fixed order).
fn quadratic_form<const N: usize>(h: &[f64; N], normal: &[[f64; N]; N]) -> f64 {
    let mut q = 0.0_f64;
    for i in 0..N {
        let mut row = 0.0_f64;
        for j in 0..N {
            row += normal[i][j] * h[j];
        }
        q += h[i] * row;
    }
    q
}

/// \\(\lVert \mathbf{v} \rVert_D = \sqrt{\sum_j (d_j v_j)^2}\\).
fn scaled_norm<const N: usize>(v: &[f64; N], d: &[f64; N]) -> f64 {
    let mut s = 0.0_f64;
    for j in 0..N {
        let t = d[j] * v[j];
        s += t * t;
    }
    s.sqrt()
}

/// Covariance at the returned point: \\(A^{-1}\\) inverted through the
/// equilibrated \\(B = D^{-1} A D^{-1}\\) for conditioning
/// (\\(A^{-1} = D^{-1} B^{-1} D^{-1}\\)). Failure is explicit — never
/// fabricated zeros.
fn covariance<const N: usize>(
    normal: &[[f64; N]; N],
    d: &[f64; N],
    d_max: f64,
) -> Result<[[f64; N]; N], CovarianceFailure> {
    let mut b = [[0.0_f64; N]; N];
    for i in 0..N {
        let di = effective_scale(d[i], d_max);
        for j in 0..N {
            b[i][j] = normal[i][j] / (di * effective_scale(d[j], d_max));
        }
    }
    let b_inv = mat_inv(&b).ok_or(CovarianceFailure::SingularNormalMatrix)?;
    let mut cov = [[0.0_f64; N]; N];
    for i in 0..N {
        let di = effective_scale(d[i], d_max);
        for j in 0..N {
            cov[i][j] = b_inv[i][j] / (di * effective_scale(d[j], d_max));
        }
    }
    let cov = mat_symmetrize(&cov);
    for i in 0..N {
        for j in 0..N {
            if !cov[i][j].is_finite() {
                return Err(CovarianceFailure::SingularNormalMatrix);
            }
        }
    }
    Ok(cov)
}

/// Nielsen rejection escalation: \\(\mu \leftarrow \max(\mu,
/// \mu_{\text{seed}}) \cdot \nu\\), \\(\nu \leftarrow 2\nu\\). Returns
/// the termination reason when the escalated \\(\mu\\) exceeds the
/// budget (the caller breaks the outer loop); `None` means retry the
/// inner loop at the new damping.
#[inline]
fn escalate_mu(mu: &mut f64, nu: &mut f64, mu_max: f64) -> Option<TerminationReason> {
    if *mu < MU_SEED {
        *mu = MU_SEED;
    }
    *mu *= *nu;
    *nu *= 2.0;
    if *mu > mu_max {
        Some(TerminationReason::DampingExhausted { mu: *mu })
    } else {
        None
    }
}

/// Construct the [`LMError::PersistentInvalidTrials`] terminal error at
/// the best (last accepted) iterate. Shared by the trial-cost and
/// rollback paths so the diagnostic payload cannot drift between them.
fn persistent_invalid_trials<E, const N: usize>(
    iteration: usize,
    consecutive: usize,
    last_source: Option<E>,
    x: &[f64; N],
    best_cost: f64,
) -> LMError<E> {
    let mut best_x = Vec::with_capacity(N);
    best_x.extend_from_slice(x);
    LMError::PersistentInvalidTrials {
        iteration,
        consecutive,
        last_source,
        best_x,
        best_cost,
    }
}

// ── Driver core ─────────────────────────────────────────────────────

/// Bookkeeping for the most recent accepted step, consumed by the
/// **cost** convergence test at the next outer iteration.
///
/// Everything here is a cross-model quantity — it only means anything
/// against the system it was measured on — which is why a
/// [`CostProblem::refresh_model`] declaration clears the whole record.
/// The step's own SIZE is not stored here for exactly that reason: it
/// stays true across a refresh, and is tracked separately so that
/// clearing this does not erase the fact that a step was accepted.
struct AcceptedStep {
    /// Whether `constrain_step` modified the step (a clamped step can
    /// never declare convergence).
    clamped: bool,
    /// Φ before the step (full-evaluation value).
    prev_cost: f64,
    /// Predicted reduction of the step.
    pred: f64,
    /// Actual reduction measured between full evaluations.
    actred: f64,
}

/// The convergence battery at an accepted point, returning the reason
/// that fired (all of which mean `converged = true`) or `None` to keep
/// iterating.
///
/// Factored out because it runs from TWO places — once per outer
/// iteration on the system the driver has been minimizing, and again
/// immediately after a [`CostProblem::refresh_model`] re-assembly,
/// which produces a different objective at the same point. Keeping one
/// implementation is what stops the two call sites from drifting into
/// different notions of "converged".
///
/// `last_accepted` gates the cost test only: pass `None` whenever the
/// last accepted step was measured against a DIFFERENT system than
/// `sys` (i.e. after a refresh), because `actred`/`pred` are then
/// cross-model quantities and comparing them is exactly what the
/// refresh exists to prevent. The gradient and step tests read only
/// `sys`, `d` and `x`, so they are always meaningful.
fn convergence_reason<const N: usize>(
    config: &LMConfig,
    x: &[f64; N],
    sys: &AssembledSystem<N>,
    d: &[f64; N],
    d_max: f64,
    last_accepted: Option<&AcceptedStep>,
) -> Option<TerminationReason> {
    // Gradient (MINPACK cosine form, multiplicative to avoid
    // division): |g_j| ≤ gtol · √A_jj · √Φ for all j, using the
    // CURRENT column norms (MINPACK semantics — the running-max d
    // would loosen the test on collapsed-sensitivity columns). At
    // an exact fit (Φ = 0) the gradient is exactly zero and the
    // test passes; a zero column forces g_j = 0, which also
    // passes. Cost is validated ≥ 0 at assembly, and the
    // is_finite guard makes a NaN threshold unreachable as a
    // matter of defense in depth — a NaN must never read as
    // "pass".
    if config.gtol > 0.0 {
        let sqrt_cost = sys.cost.sqrt();
        let mut pass = sqrt_cost.is_finite();
        for j in 0..N {
            if sys.rhs[j].abs() > config.gtol * sys.normal[j][j].sqrt() * sqrt_cost {
                pass = false;
                break;
            }
        }
        if pass {
            return Some(TerminationReason::GradientTolerance);
        }
    }

    // Step- and cost-based convergence are judged against the
    // UNDAMPED Gauss-Newton step at the accepted point. A step
    // that is tiny only because μ crushed it says nothing about
    // stationarity: on stiff valleys the rejections inflate μ by
    // orders of magnitude, and any μ-shrunken accepted step would
    // read as "converged" at an arbitrarily bad iterate (observed
    // on the 2020 CD3 capture-spanning fit at χ² ~ 1e12). The
    // undamped step is also what the legacy solver's step test
    // effectively measured (λ ≈ 1e-6 at convergence).
    //
    // A singular undamped system simply skips these tests — the
    // gradient and budget criteria still terminate.
    let h_gn = solve_damped(&sys.normal, &sys.rhs, d, d_max, 0.0)?;
    let q_gn = quadratic_form(&h_gn, &sys.normal);
    let q_pass = config.qtol > 0.0 && q_gn <= config.qtol;
    let x_pass = config.xtol > 0.0
        && scaled_norm(&h_gn, d) <= config.xtol * (scaled_norm(x, d) + config.xtol);
    if q_pass || x_pass {
        return Some(TerminationReason::StepTolerance);
    }

    // Cost tolerance (MINPACK info=1 analogue): the actual AND
    // predicted reduction of the last accepted (unclamped)
    // step both ≤ ftol·Φ_prev with a consistent model
    // (ρ ≤ 2) — AND the undamped step's own predicted gain is
    // below the same threshold, so a μ-starved step cannot
    // manufacture a plateau (for the exact GN step the
    // predicted reduction is \(\mathbf{h}^T \mathbf{g}\),
    // but the general form is used for uniformity).
    if let Some(acc) = last_accepted
        && !acc.clamped
        && config.ftol > 0.0
        && acc.prev_cost > 0.0
    {
        let threshold = config.ftol * acc.prev_cost;
        let ratio_ok = acc.pred > 0.0 && acc.actred <= 2.0 * acc.pred;
        // Same rule as the acceptance sites: when the square-root system
        // is available this must not be evaluated through A either, or the
        // convergence test is decided on a quantity with no significant
        // digits left at high row-weight dynamic range.
        let gn_pred = match (config.square_root_solve, sys.ls_rows.as_ref()) {
            (true, Some((rows, b))) => predicted_reduction_rows(&h_gn, rows, b),
            _ => predicted_reduction(&h_gn, &sys.normal, &sys.rhs),
        };
        let gn_exhausted = gn_pred <= config.ftol * sys.cost;
        if acc.actred.abs() <= threshold && acc.pred <= threshold && ratio_ok && gn_exhausted {
            return Some(TerminationReason::CostTolerance);
        }
    }

    None
}

fn solve_core<S: SystemSource<N>, const N: usize>(
    source: &mut S,
    x0: [f64; N],
    config: &LMConfig,
) -> Result<LMSolution<N>, LMError<S::Error>> {
    let mut x = x0;

    // Initial full evaluation: the only point where failure is fatal
    // (there is no accepted state to retreat to).
    let mut sys = match source.assemble(&x) {
        Ok(s) => s,
        Err(AssembleFailure::Domain(e)) => {
            return Err(LMError::InitialEvaluationFailed { source: e });
        }
        Err(AssembleFailure::NonFinite(defect)) => {
            return Err(LMError::InvalidEvaluation {
                iteration: 0,
                defect,
            });
        }
        Err(AssembleFailure::SystemNonFinite(defect)) => {
            return Err(LMError::InvalidSystem {
                iteration: 0,
                defect,
            });
        }
        Err(AssembleFailure::Hard(h)) => return Err(h.into_error(0)),
    };

    let mut d = [0.0_f64; N];
    let mut d_max = update_scaling(&mut d, &sys.normal);
    if d_max == 0.0 {
        return Err(LMError::ZeroDampingDiagonal { iteration: 0 });
    }

    // The validated initial evaluation is committed: problem-side
    // pending state (caches, diagnostics) corresponds to x0 even on
    // zero-acceptance exits, so the final-state contract holds on
    // every path.
    source.accepted(&x);

    // μ₀ = τ · max_j(A_jj / d_j²)  (≈ τ under iteration-0 scaling;
    // ≤ τ always, since d_j ≥ √A_jj — paired with the τ ≤ mu_max
    // config validation this guarantees μ₀ ≤ mu_max).
    let mut mu = {
        let mut m = 0.0_f64;
        for j in 0..N {
            let dj = effective_scale(d[j], d_max);
            let s = sys.normal[j][j] / (dj * dj);
            if s > m {
                m = s;
            }
        }
        config.tau * m
    };
    let mut nu = 2.0_f64;
    // μ as the last COMMITTED evaluation left it — the initial one, an
    // acceptance, or a re-assembly. Read only by the inexact-system path
    // (see [`CostProblem::assembly_is_exact`]): the escalation a driver
    // spends failing to improve a stand-in describes the stand-in, and is
    // discarded with it rather than inherited by the objective.
    let mut mu_committed = mu;
    // An exhaustion reached on an inexact system, carried to the top of
    // the next iteration where the forced re-assembly happens. `Some`
    // only while the system is inexact.
    let mut deferred_exhaustion: Option<TerminationReason> = None;

    let mut last_accepted: Option<AcceptedStep> = None;
    // Survives a model refresh; `last_accepted` does not (see the
    // acceptance block).
    let mut last_accepted_qnorm: Option<f64> = None;
    let mut consecutive_invalid = 0usize;
    let mut last_invalid_source: Option<S::Error> = None;

    let mut iterations = 0usize;
    let mut n_cost_evals = 0usize;
    let mut n_rejected_trials = 0usize;
    let mut n_invalid_trials = 0usize;
    let mut n_accelerated_trials = 0usize;
    let mut n_model_refreshes = 0usize;

    let mut converged = false;
    let mut reason = TerminationReason::MaxIterations;

    'outer: for iteration in 1..=config.max_iterations {
        iterations = iteration;

        // ── Convergence tests at the accepted point ──
        // Judged on the system the driver has been minimizing, with
        // the last accepted step's reductions available to the cost
        // test (they were measured against this same system).
        //
        // A criterion firing on an INEXACT system is not a verdict: what
        // it measured is the stationarity of the stand-in the problem
        // served in place of its objective, and a stand-in anchored at a
        // nearby point is stationary at its own minimizer by
        // construction. It does end the walk over that stand-in — the
        // re-assembly below is forced by it — but only the re-test on
        // what comes back may latch. See
        // [`CostProblem::assembly_is_exact`].
        let fired = convergence_reason(config, &x, &sys, &d, d_max, last_accepted.as_ref());
        let exact = source.assembly_is_exact();
        if let Some(r) = &fired
            && exact
        {
            converged = true;
            reason = r.clone();
            break 'outer;
        }

        // ── Model refresh (inexact-model problems) ──
        // A problem that evaluates a surrogate / cached linearization
        // may move its anchor HERE and nowhere else. The switch is
        // only sound together with a re-assembly: Φ(x), A, g and every
        // trial Φ(x + h) of the iteration below must describe ONE
        // objective, or ρ compares two different functions and A stops
        // being the Gauss-Newton Hessian of the compared cost. Placed
        // after the convergence tests deliberately — the tests judge
        // the objective the driver was actually minimizing, and a
        // solve that terminates here spends no refresh evaluation.
        //
        // Three triggers, one re-assembly. The problem may DECLARE the
        // refresh; or the driver forces it because the system it holds is
        // inexact and has just run out of road on that stand-in — either
        // by satisfying a convergence criterion (`fired`) or by
        // exhausting damping / inner trials at the end of the previous
        // iteration (`deferred_exhaustion`). A forced refresh is the same
        // evaluation, counted the same way; what differs is only who
        // asked for it.
        debug_assert!(
            deferred_exhaustion.is_none() || !exact,
            "an exhaustion is deferred only on an inexact system, and nothing between the \
             deferral and here assembles"
        );
        let forced = !exact && (fired.is_some() || deferred_exhaustion.is_some());
        if forced || source.refresh_model(&x) {
            // The re-assembly below is a FULL evaluation (residuals +
            // Jacobian) that no trial counter sees and that
            // `max_iterations` does not bound — for a problem whose
            // assembly is an n-body propagation it is the dominant
            // per-refresh cost. Counted here, at the declaration, so
            // the reported number is the number of assemblies the
            // refresh path actually spent (the failure arms below
            // return, so the count is only ever read on a solve that
            // completed every one of them).
            n_model_refreshes += 1;
            sys = match source.assemble(&x) {
                Ok(s) => s,
                // No retreat: the problem has already switched models,
                // so `sys` describes an objective it no longer
                // evaluates. Surface it rather than damping a stale
                // system.
                Err(AssembleFailure::Domain(e)) => {
                    return Err(LMError::ModelRefreshFailed {
                        iteration,
                        source: e,
                    });
                }
                Err(AssembleFailure::NonFinite(defect)) => {
                    return Err(LMError::InvalidEvaluation { iteration, defect });
                }
                Err(AssembleFailure::SystemNonFinite(defect)) => {
                    return Err(LMError::InvalidSystem { iteration, defect });
                }
                Err(AssembleFailure::Hard(h)) => return Err(h.into_error(iteration)),
            };
            // The whole point of a forced refresh is to obtain a system
            // the driver may judge. A second stand-in here would leave it
            // with none, and with no honest way to stop.
            if !source.assembly_is_exact() {
                return Err(LMError::InexactAssemblyPersisted { iteration });
            }
            d_max = update_scaling(&mut d, &sys.normal);
            if d_max == 0.0 {
                return Err(LMError::ZeroDampingDiagonal { iteration });
            }
            // The refreshed evaluation is the committed one at x.
            source.accepted(&x);
            // The last accepted step's reductions were measured on the
            // previous model; comparing them against this one is
            // exactly the cross-model comparison the refresh exists to
            // prevent, so the cost test restarts.
            last_accepted = None;

            // The battery above judged the model the problem just
            // replaced. The new one is a different objective at the
            // same point and may already be stationary there — a
            // surrogate refreshing onto the full model at a converged
            // iterate is precisely that case. Re-test BEFORE spending a
            // trial: otherwise the driver damps its way to
            // DampingExhausted / InnerTrialsExhausted and reports
            // failure at a point whose undamped Gauss-Newton step is
            // below `qtol`.
            //
            // Which tests re-run, and why: `last_accepted` was just
            // cleared, so the cost test is skipped BY CONSTRUCTION —
            // its actred/pred belong to the previous model and are not
            // comparable to this one. The step test (`qtol`/`xtol` on
            // the undamped Gauss-Newton step) and the MINPACK
            // gradient test read only the freshly assembled system, so
            // both are meaningful the instant the re-assembly lands.
            // The accepted DAMPED step is deliberately not a
            // convergence criterion anywhere — see the gain-ratio
            // design note; a μ-crushed step says nothing about
            // stationarity under either model.
            if let Some(r) = convergence_reason(config, &x, &sys, &d, d_max, None) {
                converged = true;
                reason = r;
                break 'outer;
            }

            // The objective is not stationary here, so an exhaustion
            // deferred from the previous iteration belonged to the
            // stand-in — and so did the damping the driver escalated
            // reaching it. Both are discarded with it: μ returns to what
            // the last committed acceptance left, and the solve carries
            // on against a system it can actually judge.
            if deferred_exhaustion.take().is_some() {
                mu = mu_committed;
                nu = 2.0;
            }
            mu_committed = mu;
        }

        // A square-root solve needs a square-root SYSTEM. There is none on
        // the normal-equations entry point ([`solve_system`]), whose caller
        // hands the driver \\(A\\) and \\(g\\) already formed — nothing was
        // squared, so there is nothing for the flag to avoid and it is
        // inert there, as documented.
        //
        // A RESIDUAL-path assembly that withheld its rows is the opposite
        // case, and absorbing it is what made this worth writing down:
        // every trial's step solve then returns `None`, the driver reads
        // each one as a rejection, and the solve ends at
        // `DampingExhausted` with μ ≈ 1e33 at its own start point, having
        // never evaluated the objective once. A step-solver outage wearing
        // a damping failure's clothes (empyrean-wv13f). Surface it.
        let square_root = config.square_root_solve && sys.ls_rows.is_some();
        if config.square_root_solve && sys.ls_rows.is_none() && sys.jacobian.is_some() {
            return Err(LMError::SquareRootSystemUnavailable { iteration });
        }

        // ── Inner trial loop: same system, escalating μ ──
        let mut accepted_this_iteration = false;
        // Set when μ escalation reached its cap. Carried out of the loop
        // rather than terminating from inside it, so that the ONE place
        // below decides what an unimproved iteration means — which
        // depends on whether the system it failed to improve was the
        // problem's objective or a stand-in for it.
        let mut exhausted: Option<TerminationReason> = None;
        for _trial in 0..config.max_inner_trials {
            // One factorization serves the velocity and (geodesic)
            // acceleration solves at this μ.
            // Square-root path: factor [J D^-1 ; sqrt(mu) I] directly.
            // Equilibration is the same D the Cholesky path uses, so the
            // two solve the identical system — one by squaring it, one
            // not. `factor` stays None here, which is what disables
            // geodesic acceleration on this path.
            let sqrt_step = if square_root {
                sys.ls_rows.as_ref().and_then(|(rows, b)| {
                    let mut qr = crate::linalg::qr::QrAccumulator::<N>::new();
                    for (row, rhs_i) in rows.iter().zip(b) {
                        let mut scaled = [0.0_f64; N];
                        for j in 0..N {
                            scaled[j] = row[j] / effective_scale(d[j], d_max);
                        }
                        qr.push_row(&scaled, *rhs_i);
                    }
                    qr.push_damping(mu);
                    let h = qr.solve().map(|h_scaled| {
                        let mut h = [0.0_f64; N];
                        for j in 0..N {
                            h[j] = h_scaled[j] / effective_scale(d[j], d_max);
                        }
                        h
                    })?;
                    // R^T R equals the SCALED damped normal matrix, so R^T
                    // is precisely the lower Cholesky factor
                    // `solve_with_factor` expects — the Givens diagonal is
                    // non-negative by construction (hypot). Handing it back
                    // keeps geodesic acceleration alive on this path with no
                    // squared matrix ever formed; without it the square-root
                    // path silently disabled the one remedy the
                    // narrow-canyon stall calls for.
                    let r = qr.r();
                    let mut l = [[0.0_f64; N]; N];
                    for (i, row) in l.iter_mut().enumerate() {
                        for (j, slot) in row.iter_mut().enumerate() {
                            *slot = r[j][i];
                        }
                    }
                    Some((l, h))
                })
            } else {
                None
            };

            let velocity = if square_root {
                sqrt_step.map(|(l, h)| (Some(l), h))
            } else {
                damped_factor(&sys.normal, &d, d_max, mu)
                    .and_then(|l| solve_with_factor(&l, &sys.rhs, &d, d_max).map(|h| (Some(l), h)))
            };
            let Some((factor, h_natural)) = velocity else {
                // Factorization failed or produced non-finite values:
                // reject (raise μ) and retry — never bail, never fall
                // back to a different solver.
                n_rejected_trials += 1;
                if let Some(r) = escalate_mu(&mut mu, &mut nu, config.mu_max) {
                    exhausted = Some(r);
                    break;
                }
                continue;
            };

            let mut h = h_natural;

            // ── Geodesic acceleration (Transtrum & Sethna 2012) ──
            // h = v + a/2 with (A + μD²) a = −Jᵀ r''_vv, where r''_vv
            // is the problem-supplied directional second derivative of
            // the residuals along v. Residual path only (needs the
            // Jacobian rows); inert when the hook returns None.
            if config.geodesic_acceleration
                && let Some(jac) = sys.jacobian.as_ref()
                && let Some(w) = source.second_directional_derivative(&x, &h_natural)
            {
                if w.len() != jac.len() {
                    return Err(LMError::DimensionMismatch {
                        iteration,
                        residuals: w.len(),
                        jacobian: jac.len(),
                    });
                }
                let mut rhs_a = [0.0_f64; N];
                for (w_i, row) in w.iter().zip(jac.iter()) {
                    for j in 0..N {
                        rhs_a[j] -= row[j] * w_i;
                    }
                }
                // Both paths supply a Cholesky factor of the damped
                // system: the Cholesky path its L, the square-root path
                // R^T, whose self-product is the scaled damped normal
                // matrix.
                if let Some(a) = factor
                    .as_ref()
                    .and_then(|l| solve_with_factor(l, &rhs_a, &d, d_max))
                {
                    let v_norm = scaled_norm(&h_natural, &d);
                    let a_norm = scaled_norm(&a, &d);
                    if std::env::var("HYPERJET_ACCEL_PROBE").is_ok() {
                        eprintln!(
                            "      ACCEL it={iteration} mu={mu:.3e} v={v_norm:.4e} a={a_norm:.4e} \
                             ratio={:.3} verdict={}",
                            a_norm / v_norm.max(1e-300),
                            if a_norm <= config.avmax * v_norm {
                                "APPLY"
                            } else {
                                "VETO"
                            },
                        );
                    }
                    if a_norm <= config.avmax * v_norm {
                        for j in 0..N {
                            h[j] += 0.5 * a[j];
                        }
                        n_accelerated_trials += 1;
                    } else {
                        // GSL avmax guard: the truncated expansion is
                        // untrustworthy here — reject through the
                        // normal μ escalation (no cost evaluation
                        // spent).
                        n_rejected_trials += 1;
                        if let Some(r) = escalate_mu(&mut mu, &mut nu, config.mu_max) {
                            exhausted = Some(r);
                            break;
                        }
                        continue;
                    }
                }
                // A failed acceleration solve simply proceeds with the
                // plain velocity step — v itself solved fine.
            }

            let h_unclamped = h;
            source.constrain(&x, &mut h);
            let clamped = h != h_unclamped;

            // Predicted reduction. For unclamped steps the model
            // decrease of the VELOCITY step is the right denominator —
            // identical to the actual step without acceleration, and
            // for accelerated steps the correction is an on-manifold
            // re-tracing the quadratic model cannot see (judging the
            // combined step against the quadratic model force-rejects
            // exactly the accelerated steps that work: on Rosenbrock
            // the combined step lands AT the minimum while the model
            // predicts an increase along it — Transtrum & Sethna 2012,
            // GSL lmaccel). A clamped step falls back to the general
            // quadratic model on the ACTUAL step, as before; pred ≤ 0
            // (or NaN) still forces rejection — no blind division
            // anywhere.
            let pred = if clamped {
                {
                    let gn = match (config.square_root_solve, sys.ls_rows.as_ref()) {
                        (true, Some((rows, b))) => predicted_reduction_rows(&h, rows, b),
                        _ => predicted_reduction(&h, &sys.normal, &sys.rhs),
                    };
                    // Fold in the exact curvature the Gauss-Newton model
                    // omits. Declines to the GN model when the problem has
                    // no second-derivative hook.
                    match (config.exact_second_order_model, sys.ls_rows.as_ref()) {
                        (true, Some((_, b))) => source
                            .second_directional_derivative(&x, &h)
                            .and_then(|w| second_order_correction(&w, b))
                            .map_or(gn, |c| gn + c),
                        _ => gn,
                    }
                }
            } else {
                match (config.square_root_solve, sys.ls_rows.as_ref()) {
                    (true, Some((rows, b))) => predicted_reduction_rows(&h_natural, rows, b),
                    _ => predicted_reduction(&h_natural, &sys.normal, &sys.rhs),
                }
            };

            let mut x_trial = [0.0_f64; N];
            for i in 0..N {
                x_trial[i] = x[i] + h[i];
            }

            n_cost_evals += 1;
            let trial_result = source.trial_cost(&x_trial);

            // Classify the trial. All comparisons are false on NaN, so
            // every Err / non-finite / model-inconsistent combination
            // lands in a rejection branch — no sentinel values.
            // NOTE: a finite trial cost does NOT reset the
            // consecutive-invalid counter yet — a provisional
            // acceptance can still be rolled back below and
            // re-classified as invalid. The counter resets only when
            // a trial is conclusively valid: a finite-cost gain-test
            // rejection, or a committed acceptance.
            let cost_trial = match trial_result {
                Ok(c) if c.is_finite() => Some(c),
                Ok(_) => {
                    n_invalid_trials += 1;
                    consecutive_invalid += 1;
                    None
                }
                Err(e) => {
                    n_invalid_trials += 1;
                    consecutive_invalid += 1;
                    last_invalid_source = Some(e);
                    None
                }
            };
            if consecutive_invalid >= config.max_consecutive_invalid {
                // Lifecycle symmetry: the trial that trips the limit
                // is still a non-committed trial — let the problem
                // discard its pending state before the error.
                source.rejected(&x_trial);
                return Err(persistent_invalid_trials(
                    iteration,
                    consecutive_invalid,
                    last_invalid_source,
                    &x,
                    sys.cost,
                ));
            }

            let accept = match cost_trial {
                Some(c) => pred > 0.0 && (sys.cost - c) > config.min_relative_decrease * pred,
                None => false,
            };

            if !accept {
                if cost_trial.is_some() {
                    n_rejected_trials += 1;
                    consecutive_invalid = 0;
                    last_invalid_source = None;
                }
                source.rejected(&x_trial);
                if let Some(r) = escalate_mu(&mut mu, &mut nu, config.mu_max) {
                    exhausted = Some(r);
                    break;
                }
                continue;
            }
            let cost_trial = cost_trial.expect("accepted trial has a finite cost");

            // Provisional acceptance: run the full evaluation at the
            // trial point BEFORE committing. A failing, non-finite, or
            // cost-inconsistent (non-decreasing — the monotonicity
            // guard) evaluation rolls the acceptance back and is
            // treated as an invalid trial.
            let assembled = source.assemble(&x_trial);
            match assembled {
                Ok(new_sys) if new_sys.cost < sys.cost => {
                    // Commit. ρ from the trial cost that drove the
                    // decision; Nielsen's smooth μ update (cube via
                    // multiplication — no libm).
                    let rho = (sys.cost - cost_trial) / pred;
                    let t = 2.0 * rho - 1.0;
                    let shrink = 1.0 - t * t * t;
                    let factor = if shrink > 1.0 / 3.0 {
                        shrink
                    } else {
                        1.0 / 3.0
                    };
                    mu *= factor;
                    // Nielsen's factor reaches 2 for marginal accepts
                    // (ρ → 0⁺), so μ can GROW through acceptances —
                    // cap it like the reject path does.
                    if mu > config.mu_max {
                        mu = config.mu_max;
                    }
                    nu = 2.0;
                    consecutive_invalid = 0;
                    last_invalid_source = None;

                    source.accepted(&x_trial);

                    // Tracked outside `last_accepted` because that
                    // record is CLEARED by a model refresh (its
                    // actred/pred stop being comparable across one),
                    // while "how big was the last step actually taken"
                    // remains true. Clearing both would make a solve
                    // that accepted many steps and then refreshed
                    // report `None` — i.e. "no step was accepted",
                    // which callers read as "the start point was
                    // already stationary".
                    last_accepted_qnorm = Some(quadratic_form(&h, &sys.normal));
                    last_accepted = Some(AcceptedStep {
                        clamped,
                        prev_cost: sys.cost,
                        pred,
                        actred: sys.cost - new_sys.cost,
                    });

                    x = x_trial;
                    sys = new_sys;
                    accepted_this_iteration = true;
                    mu_committed = mu;
                }
                Err(AssembleFailure::Hard(hard)) => {
                    // The trial never committed; discard its pending
                    // state before surfacing the contract violation.
                    source.rejected(&x_trial);
                    return Err(hard.into_error(iteration));
                }
                other => {
                    // Roll back: x unchanged, sys unchanged, no commit.
                    let source_err = match other {
                        Err(AssembleFailure::Domain(e)) => Some(e),
                        _ => None,
                    };
                    n_invalid_trials += 1;
                    consecutive_invalid += 1;
                    if let Some(e) = source_err {
                        last_invalid_source = Some(e);
                    }
                    if consecutive_invalid >= config.max_consecutive_invalid {
                        source.rejected(&x_trial);
                        return Err(persistent_invalid_trials(
                            iteration,
                            consecutive_invalid,
                            last_invalid_source,
                            &x,
                            sys.cost,
                        ));
                    }
                    source.rejected(&x_trial);
                    if let Some(r) = escalate_mu(&mut mu, &mut nu, config.mu_max) {
                        exhausted = Some(r);
                        break;
                    }
                    continue;
                }
            }

            if accepted_this_iteration {
                break;
            }
        }

        if !accepted_this_iteration {
            let stop = exhausted.unwrap_or(TerminationReason::InnerTrialsExhausted {
                trials: config.max_inner_trials,
            });
            // An INEXACT system may not end the solve either. What the
            // driver ran out of road on is the stand-in the problem
            // served, and "no step improves this linearization" is the
            // expected end of a walk across one, not a statement about
            // the objective. Defer the verdict to the top of the next
            // iteration, where the forced re-assembly turns it into a
            // question the objective can answer.
            //
            // This terminates. A system goes inexact only through an
            // assembly at an accepted point, so an exhaustion carrying a
            // stand-in had at least one acceptance behind it; and the
            // re-assembly the deferral forces is exact, so an iteration
            // that then exhausts without accepting anything stops right
            // here on the next pass.
            if source.assembly_is_exact() {
                reason = stop;
                break 'outer;
            }
            deferred_exhaustion = Some(stop);
            continue 'outer;
        }

        // Running-max scaling update at the newly accepted system.
        d_max = update_scaling(&mut d, &sys.normal);
    }

    // Every stop the driver CHOOSES already leaves `sys` exact: a
    // convergence verdict may only latch on one, and an exhaustion
    // reached on a stand-in is deferred instead of taken. The iteration
    // budget is the one stop it does not choose, so it is the one that
    // can land mid-walk — and the covariance, the gradient norm and the
    // returned q-norm below are all read off `sys`.
    if !source.assembly_is_exact() {
        n_model_refreshes += 1;
        sys = match source.assemble(&x) {
            Ok(s) => s,
            Err(AssembleFailure::Domain(e)) => {
                return Err(LMError::ModelRefreshFailed {
                    iteration: iterations,
                    source: e,
                });
            }
            Err(AssembleFailure::NonFinite(defect)) => {
                return Err(LMError::InvalidEvaluation {
                    iteration: iterations,
                    defect,
                });
            }
            Err(AssembleFailure::SystemNonFinite(defect)) => {
                return Err(LMError::InvalidSystem {
                    iteration: iterations,
                    defect,
                });
            }
            Err(AssembleFailure::Hard(h)) => return Err(h.into_error(iterations)),
        };
        if !source.assembly_is_exact() {
            return Err(LMError::InexactAssemblyPersisted {
                iteration: iterations,
            });
        }
        d_max = update_scaling(&mut d, &sys.normal);
        if d_max == 0.0 {
            return Err(LMError::ZeroDampingDiagonal {
                iteration: iterations,
            });
        }
        source.accepted(&x);
        // A deferred exhaustion was already owed this re-assembly and the
        // re-test that goes with it; the budget cut in between. Finish
        // the operation rather than reporting the budget for a walk that
        // had already ended — and rather than reporting a stall at a
        // point the objective turns out to be stationary at.
        if let Some(stop) = deferred_exhaustion.take() {
            match convergence_reason(config, &x, &sys, &d, d_max, None) {
                Some(r) => {
                    converged = true;
                    reason = r;
                }
                None => reason = stop,
            }
        }
    }

    // Final state: `sys` is the full evaluation at the returned `x` on
    // every exit path (state moved only on committed acceptances).
    let gradient_norm_scaled = {
        let mut m = 0.0_f64;
        for j in 0..N {
            let s = sys.rhs[j].abs() / effective_scale(d[j], d_max);
            if s > m {
                m = s;
            }
        }
        m
    };

    // The quantity `qtol` is tested on, evaluated at the RETURNED
    // point: the undamped Gauss-Newton step's quadratic form. Reported
    // on every exit path so a non-converged solve can be described in
    // terms of the test it failed, instead of the μ-damped accepted
    // step (which is comparable to nothing).
    let final_gn_qnorm = solve_damped(&sys.normal, &sys.rhs, &d, d_max, 0.0)
        .map(|h| quadratic_form(&h, &sys.normal));

    Ok(LMSolution {
        x,
        covariance: covariance(&sys.normal, &d, d_max),
        cost: sys.cost,
        data_cost: sys.data_cost,
        accepted_step_qnorm: last_accepted_qnorm,
        final_gn_qnorm,
        gradient_norm_scaled,
        iterations,
        n_cost_evals,
        n_rejected_trials,
        n_invalid_trials,
        n_accelerated_trials,
        n_model_refreshes,
        mu_final: mu,
        converged,
        reason,
    })
}

// ── Tests ───────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;

    /// Test domain error.
    #[derive(Debug, PartialEq)]
    struct TestError(&'static str);
    impl std::fmt::Display for TestError {
        fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
            write!(f, "{}", self.0)
        }
    }
    impl std::error::Error for TestError {}

    /// A residual problem built from a residual-vector closure, with
    /// trial costs computed from the same closure (same objective,
    /// same bits) and full lifecycle logging.
    struct Tracked<F> {
        f: F,
        accepted: Vec<Vec<f64>>,
        rejected: Vec<Vec<f64>>,
        full_evals: Vec<f64>,
    }

    impl<F> Tracked<F> {
        fn new(f: F) -> Self {
            Self {
                f,
                accepted: Vec::new(),
                rejected: Vec::new(),
                full_evals: Vec::new(),
            }
        }
    }

    impl<F, const N: usize> CostProblem<N> for Tracked<F>
    where
        F: FnMut(&[f64; N]) -> (Vec<f64>, Vec<[f64; N]>),
    {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; N]) -> Result<f64, TestError> {
            let (residuals, _) = (self.f)(x);
            Ok(residuals.iter().map(|r| r * r).sum())
        }
        fn on_step_accepted(&mut self, x: &[f64; N]) {
            self.accepted.push(x.to_vec());
        }
        fn on_step_rejected(&mut self, x_trial: &[f64; N]) {
            self.rejected.push(x_trial.to_vec());
        }
    }

    impl<F, const N: usize> ResidualProblem<N> for Tracked<F>
    where
        F: FnMut(&[f64; N]) -> (Vec<f64>, Vec<[f64; N]>),
    {
        fn evaluate(&mut self, x: &[f64; N]) -> Result<NLLSEvaluation<N>, TestError> {
            let (residuals, jacobian) = (self.f)(x);
            let cost = residuals.iter().map(|r| r * r).sum();
            self.full_evals.push(cost);
            Ok(NLLSEvaluation {
                residuals,
                jacobian,
                cost,
            })
        }
    }

    fn config() -> LMConfig {
        LMConfig::default()
    }

    /// Residual/Jacobian pair returned by 1-D test problems.
    type Resid1 = (Vec<f64>, Vec<[f64; 1]>);

    /// Rational-arithmetic overshoot problem (libm-free, so traces are
    /// bit-identical cross-platform): r(d) = d + 4d/(1+d²), d = x − 5.
    /// Around |d| ≈ 1.5 the Jacobian dips and Gauss-Newton overshoots,
    /// so the accept/reject loop must field genuine rejections.
    fn overshoot_rational(x: &[f64; 1]) -> (Vec<f64>, Vec<[f64; 1]>) {
        let d = x[0] - 5.0;
        let q = 1.0 + d * d;
        let r = d + 4.0 * d / q;
        let j = 1.0 + 4.0 * (1.0 - d * d) / (q * q);
        (vec![r], vec![[j]])
    }

    // ── Convergence on reference problems ──

    #[test]
    fn test_linear_system() {
        let mut p = Tracked::new(|x: &[f64; 2]| {
            (vec![x[0] - 3.0, x[1] - 7.0], vec![[1.0, 0.0], [0.0, 1.0]])
        });
        let mut cfg = config();
        cfg.xtol = 1e-14;
        let sol = solve(&mut p, [0.0; 2], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 3.0).abs() < 1e-10, "x0={}", sol.x[0]);
        assert!((sol.x[1] - 7.0).abs() < 1e-10, "x1={}", sol.x[1]);
        let cov = sol.covariance.unwrap();
        assert!((cov[0][0] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_overdetermined_linear() {
        let xs = [0.0, 1.0, 2.0, 3.0];
        let ys = [1.0, 3.0, 5.0, 7.0];
        let mut p = Tracked::new(move |p: &[f64; 2]| {
            let residuals: Vec<f64> = xs
                .iter()
                .zip(ys.iter())
                .map(|(&x, &y)| p[0] * x + p[1] - y)
                .collect();
            let jacobian: Vec<[f64; 2]> = xs.iter().map(|&x| [x, 1.0]).collect();
            (residuals, jacobian)
        });
        let sol = solve(&mut p, [0.0; 2], &config(), None).unwrap();
        assert!(sol.converged);
        assert!((sol.x[0] - 2.0).abs() < 1e-8, "a={}", sol.x[0]);
        assert!((sol.x[1] - 1.0).abs() < 1e-8, "b={}", sol.x[1]);
    }

    fn rosenbrock(x: &[f64; 2]) -> (Vec<f64>, Vec<[f64; 2]>) {
        (
            vec![10.0 * (x[1] - x[0] * x[0]), 1.0 - x[0]],
            vec![[-20.0 * x[0], 10.0], [-1.0, 0.0]],
        )
    }

    #[test]
    fn test_rosenbrock() {
        let mut p = Tracked::new(rosenbrock);
        let sol = solve(&mut p, [-1.0, 1.0], &config(), None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 1.0).abs() < 1e-6, "x0={}", sol.x[0]);
        assert!((sol.x[1] - 1.0).abs() < 1e-6, "x1={}", sol.x[1]);
    }

    #[test]
    fn test_circle_fit() {
        let angles: [f64; 8] = [0.0, 0.7, 1.4, 2.1, 2.8, 3.5, 4.2, 4.9];
        let data: Vec<(f64, f64)> = angles
            .iter()
            .map(|&a| (2.0 + 5.0 * a.cos(), 3.0 + 5.0 * a.sin()))
            .collect();
        let mut p = Tracked::new(move |p: &[f64; 3]| {
            let (cx, cy, r) = (p[0], p[1], p[2]);
            let mut residuals = Vec::new();
            let mut jacobian = Vec::new();
            for &(x, y) in &data {
                let dx = x - cx;
                let dy = y - cy;
                let dist = (dx * dx + dy * dy).sqrt();
                residuals.push(dist - r);
                jacobian.push([-dx / dist, -dy / dist, -1.0]);
            }
            (residuals, jacobian)
        });
        let sol = solve(&mut p, [0.0, 0.0, 1.0], &config(), None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 2.0).abs() < 1e-7, "cx={}", sol.x[0]);
        assert!((sol.x[1] - 3.0).abs() < 1e-7, "cy={}", sol.x[1]);
        assert!((sol.x[2] - 5.0).abs() < 1e-7, "r={}", sol.x[2]);
    }

    /// The motivating defect class (empyrean-ju91): far from the
    /// minimum of r = ln(1 + (x−5)²) the Jacobian collapses and pure
    /// Gauss-Newton overshoots by orders of magnitude. The removed
    /// always-accept solver provably diverged here without a
    /// problem-supplied step clamp; the accept/reject loop must
    /// converge with NO clamp.
    #[test]
    fn test_overshoot_converges_without_clamp() {
        let mut p = Tracked::new(|x: &[f64; 1]| {
            let d = x[0] - 5.0;
            let r = (1.0 + d * d).ln();
            let j = 2.0 * d / (1.0 + d * d);
            (vec![r], vec![[j]])
        });
        let mut cfg = config();
        cfg.max_iterations = 300;
        let sol = solve(&mut p, [100.0], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 5.0).abs() < 1e-3, "x={}", sol.x[0]);
        assert!(
            sol.n_rejected_trials > 0,
            "overshoot should exercise rejections"
        );
    }

    // ── τ-insensitivity: the dead-zone regression ──

    /// A correct accept/reject solver is insensitive to the initial
    /// damping choice. This encodes the ju91/CD3 dead-zone symptom
    /// (convergence non-monotonic in λ_initial) as a permanent
    /// invariant.
    #[test]
    fn test_tau_insensitivity_sweep() {
        let taus = [1e-8, 1e-6, 1e-4, 1e-3, 1e-1, 1.0, 1e2];
        for &tau in &taus {
            let mut cfg = config();
            cfg.tau = tau;
            cfg.max_iterations = 500;

            let mut p = Tracked::new(rosenbrock);
            let sol = solve(&mut p, [-1.0, 1.0], &cfg, None).unwrap();
            assert!(sol.converged, "rosenbrock tau={tau}: {:?}", sol.reason);
            assert!(
                (sol.x[0] - 1.0).abs() < 1e-6 && (sol.x[1] - 1.0).abs() < 1e-6,
                "rosenbrock tau={tau}: x={:?}",
                sol.x
            );

            let mut p = Tracked::new(overshoot_rational);
            let sol = solve(&mut p, [6.5], &cfg, None).unwrap();
            assert!(sol.converged, "overshoot tau={tau}: {:?}", sol.reason);
            assert!(
                (sol.x[0] - 5.0).abs() < 1e-4,
                "overshoot tau={tau}: x={}",
                sol.x[0]
            );

            let mut p = Tracked::new(|x: &[f64; 1]| {
                let d = x[0] - 5.0;
                let r = (1.0 + d * d).ln();
                let j = 2.0 * d / (1.0 + d * d);
                (vec![r], vec![[j]])
            });
            let sol = solve(&mut p, [100.0], &cfg, None).unwrap();
            assert!(sol.converged, "log-overshoot tau={tau}: {:?}", sol.reason);
            assert!(
                (sol.x[0] - 5.0).abs() < 1e-3,
                "log-overshoot tau={tau}: x={}",
                sol.x[0]
            );
        }
    }

    // ── Prior handling ──

    #[test]
    fn test_prior_map_solution() {
        // r = x − 10 with prior N(0, 1): MAP at x = 5, covariance 1/2.
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0] - 10.0], vec![[1.0]]));
        let prior = NLLSPrior {
            mean: [0.0],
            covariance_inv: [[1.0]],
        };
        let sol = solve(&mut p, [0.0], &config(), Some(&prior)).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 5.0).abs() < 1e-8, "x={}", sol.x[0]);
        let cov = sol.covariance.unwrap();
        assert!((cov[0][0] - 0.5).abs() < 1e-10, "cov={}", cov[0][0]);
        // cost = data + prior penalty; data_cost = data only.
        assert!((sol.cost - 50.0).abs() < 1e-6, "cost={}", sol.cost);
        assert!(
            (sol.data_cost - 25.0).abs() < 1e-6,
            "data={}",
            sol.data_cost
        );
    }

    /// Defect-13 regression: with a weak data term and a tight prior,
    /// the solver must land on the MAP point, not the data minimum —
    /// the prior penalty enters both compared costs, so steps that
    /// trade a small data gain for a large prior penalty are rejected.
    #[test]
    fn test_prior_penalty_in_acceptance() {
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![0.1 * (x[0] - 10.0)], vec![[0.1]]));
        let prior = NLLSPrior {
            mean: [0.0],
            covariance_inv: [[1.0]],
        };
        let mut cfg = config();
        cfg.xtol = 1e-13;
        // Start AT the data minimum: every step toward the MAP point
        // INCREASES the data cost, so a data-only acceptance objective
        // (the defect) would reject every trial and never leave
        // x = 10 — this start point is what makes the test
        // discriminating.
        let sol = solve(&mut p, [10.0], &cfg, Some(&prior)).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        // MAP: minimize 0.01(x−10)² + x² → x = 0.1/1.01.
        let expected = 0.1 / 1.01;
        assert!(
            (sol.x[0] - expected).abs() < 1e-6,
            "x={} expected={expected}",
            sol.x[0]
        );
    }

    /// Uninformative data (zero Jacobian, finite residuals) with a
    /// prior: the posterior IS the prior. The solver must converge to
    /// the prior mean with the prior covariance — honestly, with
    /// finite cost (the old solver's INF-sentinel false convergence is
    /// structurally impossible here).
    #[test]
    fn test_prior_zero_jacobian_posterior_is_prior() {
        let mut p = Tracked::new(|_x: &[f64; 1]| (vec![5.0], vec![[0.0]]));
        let prior = NLLSPrior {
            mean: [2.0],
            covariance_inv: [[1.0]],
        };
        let sol = solve(&mut p, [0.0], &config(), Some(&prior)).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 2.0).abs() < 1e-8, "x={}", sol.x[0]);
        let cov = sol.covariance.unwrap();
        assert!((cov[0][0] - 1.0).abs() < 1e-8, "cov={}", cov[0][0]);
        assert!(sol.cost.is_finite());
        assert!((sol.data_cost - 25.0).abs() < 1e-10);
    }

    // ── Invalid inputs and error axes ──

    #[test]
    fn test_initial_domain_error() {
        struct Failing;
        impl CostProblem<1> for Failing {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Err(TestError("cost"))
            }
        }
        impl ResidualProblem<1> for Failing {
            fn evaluate(&mut self, _x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Err(TestError("propagation failed"))
            }
        }
        let err = solve(&mut Failing, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(err, LMError::InitialEvaluationFailed { .. }),
            "{err:?}"
        );
    }

    #[test]
    fn test_nonfinite_initial_cost_is_error() {
        // An INF cost at x0 must be a loud error, not a sentinel the
        // solver silently converges on.
        struct InfCost;
        impl CostProblem<1> for InfCost {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(f64::INFINITY)
            }
        }
        impl ResidualProblem<1> for InfCost {
            fn evaluate(&mut self, _x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![0.0],
                    jacobian: vec![[0.0]],
                    cost: f64::INFINITY,
                })
            }
        }
        let err = solve(&mut InfCost, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidEvaluation {
                    iteration: 0,
                    defect: EvaluationDefect::Cost
                }
            ),
            "{err:?}"
        );
    }

    #[test]
    fn test_empty_residuals_and_dimension_mismatch() {
        let mut p = Tracked::new(|_x: &[f64; 1]| (vec![], vec![]));
        let err = solve(&mut p, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(err, LMError::EmptyResiduals { iteration: 0 }),
            "{err:?}"
        );

        let mut p = Tracked::new(|_x: &[f64; 1]| (vec![1.0, 2.0], vec![[1.0]]));
        let err = solve(&mut p, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::DimensionMismatch {
                    iteration: 0,
                    residuals: 2,
                    jacobian: 1
                }
            ),
            "{err:?}"
        );
    }

    /// The exact second-order correction must reproduce the TRUE model
    /// reduction on a large-residual problem, where the Gauss-Newton model
    /// does not.
    ///
    /// Residuals r_i(x) = c_i - x^2 with c_i chosen so the fit cannot reach
    /// zero: at the solution the residuals stay large, which is exactly the
    /// regime where S = sum_i r_i grad^2 r_i stops being negligible. Here
    /// grad^2 r_i = -2 for every i, so S = -2 * sum_i r_i is known in closed
    /// form and the correction can be checked against arithmetic rather than
    /// against another implementation.
    #[test]
    fn test_second_order_correction_matches_the_closed_form_s_term() {
        // r_i = c_i - x^2, J_i = -2x, d2r_i/dx2 = -2.
        let c = [1.0_f64, 4.0, 9.0];
        let x = [1.5_f64];
        let h = [0.25_f64];

        let r: Vec<f64> = c.iter().map(|ci| ci - x[0] * x[0]).collect();
        let b: Vec<f64> = r.iter().map(|ri| -ri).collect();
        // Second directional derivative along h: h^T grad^2 r_i h = -2 h^2.
        let w: Vec<f64> = c.iter().map(|_| -2.0 * h[0] * h[0]).collect();

        let correction = second_order_correction(&w, &b).expect("lengths match");

        // h^T S h = sum_i r_i * (h^T grad^2 r_i h), and the corrected model
        // reduction is pred - h^T S h = pred + correction.
        let hsh: f64 = r.iter().zip(&w).map(|(ri, wi)| ri * wi).sum();
        assert!(
            (correction + hsh).abs() < 1e-12,
            "correction {correction} should be -(h^T S h) = {}",
            -hsh
        );

        // And S is not negligible NEXT TO the Gauss-Newton term, which is the
        // statement that matters: an absolute threshold would only measure the
        // fixture's units. J_i = -2x for every row, so h^T J^T J h is exact
        // arithmetic here too. If this ratio ever collapses, the fixture has
        // stopped exercising the large-residual regime and the test is vacuous.
        let jtj: f64 = c.iter().map(|_| 4.0 * x[0] * x[0]).sum();
        let hjjh = jtj * h[0] * h[0];
        assert!(
            hsh.abs() / hjjh > 0.1,
            "S contributes only {:.1}% of the Gauss-Newton term; the fixture is \
             no longer large-residual",
            100.0 * hsh.abs() / hjjh,
        );
    }

    /// A length mismatch yields no correction rather than a partial one.
    /// Applying curvature from some residuals and not others would be a
    /// quietly wrong model, which is worse than the Gauss-Newton model it
    /// was meant to improve.
    #[test]
    fn test_second_order_correction_refuses_a_partial_application() {
        assert!(second_order_correction(&[1.0, 2.0], &[1.0]).is_none());
    }

    /// THE CASE THE SQUARE-ROOT PATH EXISTS FOR.
    ///
    /// A linear least-squares problem whose rows carry weights differing by
    /// 1e7. The problem itself is well posed and exactly solvable — Van Loan
    /// 1985 — but the normal matrix has condition number ~1e14 to 1e16, so
    /// the Cholesky path either fails to converge or lands far from the
    /// answer, while the QR path never squares anything and gets it.
    ///
    /// This is the shape of a joint optical + radar orbit fit, where
    /// pre-weighted radar delay rows carry |J|/sigma five orders above
    /// optical's.
    #[test]
    fn test_square_root_solve_survives_stiff_row_weights() {
        const STIFF: f64 = 1.0e7;
        // Exact solution [1, -2]. Two stiff rows and two ordinary ones,
        // consistent, so the residual minimum is zero.
        struct Stiff;
        impl CostProblem<2> for Stiff {
            type Error = std::convert::Infallible;
            fn evaluate_cost(&mut self, x: &[f64; 2]) -> Result<f64, Self::Error> {
                let r = Self::residuals(x);
                Ok(r.iter().map(|v| v * v).sum())
            }
        }
        impl Stiff {
            fn residuals(x: &[f64; 2]) -> [f64; 4] {
                [
                    STIFF * (x[0] + x[1] - (-1.0)),
                    STIFF * (x[0] - x[1] - 3.0),
                    x[0] - 1.0,
                    x[1] - (-2.0),
                ]
            }
        }
        impl ResidualProblem<2> for Stiff {
            fn evaluate(&mut self, x: &[f64; 2]) -> Result<NLLSEvaluation<2>, Self::Error> {
                let r = Self::residuals(x);
                let jac = vec![[STIFF, STIFF], [STIFF, -STIFF], [1.0, 0.0], [0.0, 1.0]];
                Ok(NLLSEvaluation {
                    residuals: r.to_vec(),
                    jacobian: jac,
                    cost: r.iter().map(|v| v * v).sum(),
                })
            }
        }

        let solve_with = |square_root: bool| {
            let config = LMConfig {
                square_root_solve: square_root,
                max_iterations: 100,
                ..Default::default()
            };
            solve(&mut Stiff, [0.0, 0.0], &config, None)
        };

        let sqrt_sol = solve_with(true).expect("square-root path must solve");
        let err_sqrt = (sqrt_sol.x[0] - 1.0).abs().max((sqrt_sol.x[1] + 2.0).abs());
        assert!(
            err_sqrt < 1e-6,
            "square-root path landed at {:?}, want [1, -2] (err {err_sqrt:.3e})",
            sqrt_sol.x,
        );

        // The control: the same problem through the normal equations. It is
        // allowed to fail outright; what it must NOT do is quietly beat the
        // square-root path, which would make this test vacuous.
        let cholesky_err = solve_with(false)
            .ok()
            .map(|s| (s.x[0] - 1.0).abs().max((s.x[1] + 2.0).abs()));
        match cholesky_err {
            None => { /* failed outright — the motivating symptom */ }
            Some(e) => assert!(
                e >= err_sqrt,
                "normal-equations path was MORE accurate ({e:.3e}) than the \
                 square-root path ({err_sqrt:.3e}); this fixture no longer \
                 exercises the conditioning difference it was written for"
            ),
        }
    }

    /// With the flag off, nothing moves. The square-root path is opt-in and
    /// must leave every existing fit bit-identical.
    #[test]
    fn test_square_root_solve_is_inert_when_disabled() {
        let off = LMConfig::default();
        assert!(!off.square_root_solve, "must default off");
    }

    /// Streaming 19,040 rows through the Givens QR must not accumulate
    /// error, and the damped factorization must stay a true factorization
    /// all the way to \\(\mu = 10^{35}\\).
    ///
    /// # Why this size and why this shape
    ///
    /// The square-root path was gated on 3-row systems. When a production
    /// orbit fit then stalled — 9,520 optical observations over a
    /// 6,276-day arc, two rows each, damping escalated to
    /// \\(\mu \approx 10^{35}\\), `DampingExhausted` — the leading
    /// hypothesis was accumulated rotation error over those 19,040 rows
    /// (empyrean-wv13f). It was not: measured on the real assembled
    /// system, \\(R^\top R\\) reproduced the scaled damped normal matrix
    /// to 1.4e-14 relative at every \\(\mu\\) from 0 to \\(10^{40}\\), and
    /// the two step solves agreed to 1e-13. This pins that, so the
    /// hypothesis cannot come back.
    ///
    /// The system is consistent with a KNOWN exact solution, so accuracy
    /// is measured against arithmetic rather than against the other
    /// solver.
    ///
    /// # What makes a stiff fixture stiff
    ///
    /// Worth stating, because two earlier drafts of this fixture were
    /// vacuous. Neither column-scale spread nor row-weight spread survives
    /// the Moré equilibration \\(B = D^{-1} A D^{-1}\\) that BOTH paths
    /// already apply: eight decades of column units gave
    /// \\(\kappa(A) = 10^{16}\\) and \\(\kappa(B) = 1.05\\), and the
    /// Cholesky path solved it more accurately than the QR. Equilibration
    /// normalizes each column by its own scale, so any stiffness that is
    /// only a choice of units is gone before either solver runs.
    ///
    /// What equilibration cannot remove is near-COLLINEARITY: two
    /// parameters the data separates only through a part-per-million
    /// difference. That is the shape an orbit fit really has — the
    /// along-track and non-gravitational directions are distinguished by a
    /// small residue on top of a large common signal — and it is the case
    /// where squaring genuinely destroys the answer. Here it leaves
    /// \\(\kappa(B) = 4.0\times 10^{12}\\): the QR returns the exact
    /// solution to 1.1e-10 while the Cholesky path reaches only 4.3e-3,
    /// seven orders apart.
    #[test]
    fn test_square_root_solve_streams_a_production_scale_system() {
        const M: usize = 19_040;
        const N: usize = 6;
        // Two nearly parallel columns: the data separates the last two
        // parameters only through a part-per-SEPARATION difference, so the
        // information distinguishing them is a small residue on top of a
        // large common component. Squaring squares that residue.
        const SEPARATION: f64 = 1e-6;
        let exact = [2.0_f64, -3.0, 0.5, 7.0, -1.25, 4.0];

        // Deterministic, dependency-free row generator (xorshift64).
        let mut state = 0x2545_F491_4F6C_DD1D_u64;
        let mut uniform = move || {
            state ^= state << 13;
            state ^= state >> 7;
            state ^= state << 17;
            (state >> 11) as f64 / (1_u64 << 53) as f64 - 0.5
        };
        let mut rows = Vec::with_capacity(M);
        let mut rhs = Vec::with_capacity(M);
        for _ in 0..M {
            let mut row = [0.0_f64; N];
            for slot in row.iter_mut() {
                *slot = uniform();
            }
            row[N - 1] = row[N - 2] + SEPARATION * uniform();
            let b: f64 = (0..N).map(|j| row[j] * exact[j]).sum();
            rows.push(row);
            rhs.push(b);
        }

        // The undamped solve must return the exact solution. Column
        // equilibration is what the driver applies, so apply it here too.
        let mut normal = [[0.0_f64; N]; N];
        let mut grad = [0.0_f64; N];
        for (row, b_i) in rows.iter().zip(&rhs) {
            for j in 0..N {
                for k in 0..N {
                    normal[j][k] += row[j] * row[k];
                }
                grad[j] += row[j] * b_i;
            }
        }
        let normal = mat_symmetrize(&normal);
        let mut d = [0.0_f64; N];
        let d_max = update_scaling(&mut d, &normal);

        use crate::linalg::qr::QrAccumulator;
        let factored = |mu: f64| -> Option<(QrAccumulator<N>, [f64; N])> {
            let mut qr = QrAccumulator::<N>::new();
            for (row, b_i) in rows.iter().zip(&rhs) {
                let mut scaled = [0.0_f64; N];
                for j in 0..N {
                    scaled[j] = row[j] / effective_scale(d[j], d_max);
                }
                qr.push_row(&scaled, *b_i);
            }
            qr.push_damping(mu);
            let h = qr.solve().map(|hs| {
                let mut h = [0.0_f64; N];
                for j in 0..N {
                    h[j] = hs[j] / effective_scale(d[j], d_max);
                }
                h
            })?;
            Some((qr, h))
        };

        // State the conditioning the gate actually runs at, AFTER the
        // equilibration both paths apply — the only conditioning number
        // that means anything here. Without this a future edit could
        // quietly turn the fixture benign and every tolerance below would
        // stop measuring anything.
        let mut equilibrated = [[0.0_f64; N]; N];
        for i in 0..N {
            for j in 0..N {
                equilibrated[i][j] =
                    normal[i][j] / (effective_scale(d[i], d_max) * effective_scale(d[j], d_max));
            }
        }
        let (eigenvalues, _) = mat_symmetric_eigen(&equilibrated).expect("symmetric");
        let lo = eigenvalues.iter().copied().fold(f64::INFINITY, f64::min);
        let hi = eigenvalues.iter().copied().fold(0.0_f64, f64::max);
        let kappa_b = hi / lo;
        assert!(
            kappa_b > 1e11,
            "fixture drifted: kappa(D^-1 A D^-1) = {kappa_b:.3e}. Equilibration has \
             absorbed the weighting, so forming the normal equations no longer costs \
             anything and this gate is vacuous"
        );

        let (_, h) = factored(0.0).expect("the undamped system is full rank");
        let mut worst_sqrt = 0.0_f64;
        for j in 0..N {
            let err = (h[j] - exact[j]).abs() / exact[j].abs();
            worst_sqrt = worst_sqrt.max(err);
            assert!(
                err < 1e-8,
                "component {j} after {M} streamed rows: {} vs exact {} \
                 (relative error {err:.3e}) — the Givens accumulation has drifted",
                h[j],
                exact[j],
            );
        }

        // The control, and the reason the whole path exists: the same
        // system through the normal equations. It may fail outright; what
        // it must NOT do is match the square-root path, because then
        // kappa(A) has drifted somewhere f64 can still cope with and the
        // gate has stopped measuring anything.
        let cholesky_worst = solve_damped(&normal, &grad, &d, d_max, 0.0).map(|h| {
            (0..N).fold(0.0_f64, |m, j| {
                m.max((h[j] - exact[j]).abs() / exact[j].abs())
            })
        });
        match cholesky_worst {
            None => { /* the factorization gave up — the motivating symptom */ }
            Some(e) => assert!(
                e > 1e3 * worst_sqrt,
                "normal equations reached {e:.3e} against the square-root path's \
                 {worst_sqrt:.3e}; at kappa(B) = {kappa_b:.3e} they should be \
                 orders apart, so this fixture is no longer stiff"
            ),
        }

        // R^T R must remain a factorization of the SCALED damped normal
        // matrix at every damping, including where mu*I swamps the data
        // entirely.
        for mu in [0.0, 1.0, 1e8, 1e20, 1e35] {
            let (qr, h) = factored(mu).expect("a damped system is always full rank");
            assert!(
                scaled_norm(&h, &d).is_finite(),
                "mu {mu:.0e} produced a non-finite step {h:?}"
            );
            let r = qr.r();
            for i in 0..N {
                for j in 0..N {
                    let rtr: f64 = (0..N).map(|k| r[k][i] * r[k][j]).sum();
                    let want = normal[i][j]
                        / (effective_scale(d[i], d_max) * effective_scale(d[j], d_max))
                        + if i == j { mu } else { 0.0 };
                    let err = (rtr - want).abs() / want.abs().max(1e-300);
                    assert!(
                        err < 1e-12,
                        "mu {mu:.0e}: (R^T R)[{i}][{j}] = {rtr} != {want} \
                         (relative {err:.3e})"
                    );
                }
            }
        }

        // And the predicted reduction taken from the rows agrees with the
        // one taken from A and g — the quantity the gain ratio divides by.
        let (_, h) = factored(1.0).expect("damped");
        let from_rows = predicted_reduction_rows(&h, &rows, &rhs);
        let from_normal = predicted_reduction(&h, &normal, &grad);
        assert!(
            (from_rows - from_normal).abs() <= 1e-8 * from_normal.abs(),
            "predicted reduction disagrees between the rows ({from_rows}) and \
             the normal equations ({from_normal})"
        );
    }

    /// A linear least-squares problem with a deliberately partial prior:
    /// information on \\(x_0\\), none on \\(x_1\\).
    struct PartialPrior;

    impl PartialPrior {
        fn residuals(x: &[f64; 2]) -> [f64; 3] {
            [x[0] - 1.0, x[1] - 2.0, 0.5 * x[0] + 0.5 * x[1] - 1.5]
        }
        fn jacobian() -> Vec<[f64; 2]> {
            vec![[1.0, 0.0], [0.0, 1.0], [0.5, 0.5]]
        }
    }

    impl CostProblem<2> for PartialPrior {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 2]) -> Result<f64, TestError> {
            Ok(Self::residuals(x).iter().map(|v| v * v).sum())
        }
    }

    impl ResidualProblem<2> for PartialPrior {
        fn evaluate(&mut self, x: &[f64; 2]) -> Result<NLLSEvaluation<2>, TestError> {
            let r = Self::residuals(x);
            Ok(NLLSEvaluation {
                cost: r.iter().map(|v| v * v).sum(),
                residuals: r.to_vec(),
                jacobian: Self::jacobian(),
            })
        }
    }

    /// A PARTIAL prior — precision on one component, none on the other —
    /// must solve the same on both paths.
    ///
    /// Such a precision matrix is positive SEMI-definite, and semi-definite
    /// matrices have no Cholesky factor. The normal-equations path never
    /// cared, because it only ever ADDS \\(P_0^{-1}\\) to \\(A\\); the
    /// square-root path needs a factor, and a Cholesky-only construction
    /// withheld the whole square-root system instead. The failure that
    /// produced was not a visibly missing prior — it was every step solve
    /// returning `None`, read as a trial rejection, μ escalating to
    /// `DampingExhausted` at 1.3e33, and a solve reporting a damping
    /// failure at its own start point with ZERO cost evaluations spent
    /// (empyrean-wv13f).
    ///
    /// The assertion is therefore on agreement AND on the solve having
    /// done work: an identical-to-the-Cholesky-path answer is the
    /// statement, and `n_cost_evals > 0` is what distinguishes solving the
    /// problem from never having started.
    #[test]
    fn test_square_root_solve_accepts_a_rank_deficient_prior() {
        // Information along ONE direction and none across it, and that
        // direction is deliberately OFF-AXIS: P = w wᵀ for a unit w at
        // 30°. `validate_prior` admits it (only a NEGATIVE diagonal is a
        // defect) — a zero precision means "no information", not
        // "invalid".
        //
        // The off-axis part is load-bearing. An axis-aligned prior like
        // [[1,0],[0,0]] is symmetric under transposing its own
        // eigenvector matrix, so it cannot distinguish rows built from
        // \\(V\\)'s columns from rows built from its rows — the exact
        // orientation error that would leave \\(S^\top S = V^\top \Lambda
        // V\\) instead of \\(V \Lambda V^\top\\), producing finite,
        // plausible, WRONG cross-terms. Anisotropy is what makes this
        // test able to fail.
        let (c, s) = (0.75_f64.sqrt(), 0.5_f64);
        let prior = NLLSPrior {
            mean: [0.0, 0.0],
            covariance_inv: [[c * c, c * s], [c * s, s * s]],
        };

        let solve_with = |square_root: bool| {
            let cfg = LMConfig {
                square_root_solve: square_root,
                ..config()
            };
            let mut p = PartialPrior;
            solve(&mut p, [5.0, -5.0], &cfg, Some(&prior)).expect("both paths must solve")
        };

        let ne = solve_with(false);
        let sq = solve_with(true);

        assert!(
            sq.converged,
            "square-root path did not converge on a semi-definite prior: {:?}",
            sq.reason
        );
        assert!(
            sq.n_cost_evals > 0,
            "square-root path terminated without evaluating the objective once — \
             the step solve is failing on every trial and μ escalation is wearing \
             the blame"
        );
        for i in 0..2 {
            assert!(
                (sq.x[i] - ne.x[i]).abs() <= 1e-9 * ne.x[i].abs().max(1.0),
                "component {i}: square-root {} vs normal equations {}",
                sq.x[i],
                ne.x[i]
            );
        }
        // The prior is genuinely in the answer: without it the fit would
        // sit at the unpenalized least-squares solution.
        let mut unprior = PartialPrior;
        let free = solve(&mut unprior, [5.0, -5.0], &config(), None).expect("unpriored solve");
        assert!(
            (sq.x[0] - free.x[0]).abs() > 1e-3,
            "the prior left no trace on x0 ({} vs unpriored {}); the fixture has \
             stopped testing that the prior rows are present at all",
            sq.x[0],
            free.x[0],
        );
    }

    /// [`prior_rows`] returns a genuine square root: \\(S^\top S =
    /// P_0^{-1}\\) exactly, definite or not, and in a basis the precision
    /// matrix does not share with the coordinate axes.
    #[test]
    fn test_prior_rows_reproduce_the_precision_matrix() {
        // Rank 2 of 3, off-axis: P = 4 u u^T + 1 v v^T with u, v
        // orthonormal and neither an axis. The third direction carries no
        // information at all.
        let u = [0.6, 0.8, 0.0];
        let v = [0.0, 0.0, 1.0];
        let mut p = [[0.0_f64; 3]; 3];
        for i in 0..3 {
            for j in 0..3 {
                p[i][j] = 4.0 * u[i] * u[j] + v[i] * v[j];
            }
        }
        let rows = prior_rows(&p).expect("a semi-definite precision has a square root");
        assert!(
            rows.len() <= 3,
            "at most one row per direction, got {rows:?}"
        );
        for i in 0..3 {
            for j in 0..3 {
                let sts: f64 = rows.iter().map(|r| r[i] * r[j]).sum();
                assert!(
                    (sts - p[i][j]).abs() < 1e-12,
                    "S^T S [{i}][{j}] = {sts} != {}",
                    p[i][j]
                );
            }
        }
        // The unconstrained direction carries no penalty: w = u × v is
        // orthogonal to both, so ‖S w‖ must be at the eigensolver's noise
        // floor and not at the scale of the prior itself.
        let w = [
            u[1] * v[2] - u[2] * v[1],
            u[2] * v[0] - u[0] * v[2],
            u[0] * v[1] - u[1] * v[0],
        ];
        let sw: f64 = rows
            .iter()
            .map(|r| {
                let d: f64 = (0..3).map(|j| r[j] * w[j]).sum();
                d * d
            })
            .sum::<f64>()
            .sqrt();
        assert!(
            sw < 1e-7,
            "the unconstrained direction picked up a penalty ‖Sw‖ = {sw:.3e}"
        );

        // A precision that is already diagonal has an exact zero
        // eigenvalue, and an exact zero contributes no row at all.
        let axis = [[1.0_f64, 0.0], [0.0, 0.0]];
        assert_eq!(
            prior_rows(&axis).expect("semi-definite").len(),
            1,
            "an exactly-zero precision direction must contribute no row"
        );

        // A positive-definite precision still goes through the Cholesky,
        // whose rows are exactly L^T — the pre-existing construction, kept
        // bit-for-bit so no definite prior's fit moves.
        let pd = [[4.0, 1.0], [1.0, 3.0]];
        let l = mat_cholesky(&pd).expect("definite");
        let rows = prior_rows(&pd).expect("definite");
        assert_eq!(rows.len(), 2);
        for i in 0..2 {
            for j in 0..2 {
                assert_eq!(rows[i][j], l[j][i], "row {i} col {j} is not L^T");
            }
        }
    }

    /// A precision matrix with a genuinely NEGATIVE eigenvalue is not a
    /// precision matrix, and no set of rows reproduces its penalty — it is
    /// not a sum of squares in any basis. The square-root path must SAY so.
    ///
    /// The behaviour being pinned is the diagnosis, not the refusal. Left
    /// unreported, the same condition presents as
    /// `DampingExhausted` at astronomical μ, which sends the reader to the
    /// damping policy instead of to the prior.
    #[test]
    fn test_square_root_solve_reports_a_prior_it_cannot_factor() {
        // Indefinite: eigenvalues 3 and -1, both diagonal entries
        // positive, so `validate_prior` (which reads only the diagonal)
        // admits it.
        let prior = NLLSPrior {
            mean: [0.0, 0.0],
            covariance_inv: [[1.0, 2.0], [2.0, 1.0]],
        };
        let cfg = LMConfig {
            square_root_solve: true,
            ..config()
        };
        let mut p = PartialPrior;
        let err = solve(&mut p, [5.0, -5.0], &cfg, Some(&prior))
            .expect_err("an indefinite prior has no square-root form");
        assert!(
            matches!(err, LMError::SquareRootSystemUnavailable { .. }),
            "expected the withheld square-root system to be named, got {err:?}"
        );
    }

    /// `square_root_solve` on the normal-equations entry point is INERT,
    /// as documented — not fatal.
    ///
    /// [`solve_system`] callers hand the driver \\(A\\) and \\(g\\)
    /// directly, so there is no squaring left for a square-root solve to
    /// avoid and no least-squares rows to factor. A consumer that sets the
    /// flag once in a shared config and reaches this entry point must get
    /// its fit, not a damping failure at its start point.
    #[test]
    fn test_square_root_solve_is_inert_on_the_system_path() {
        let cfg = LMConfig {
            square_root_solve: true,
            ..config()
        };
        let mut p = QuadSystem {
            q: [[2.0, 0.0], [0.0, 8.0]],
            a: [1.0, 2.0],
            c: 3.0,
        };
        let sol = solve_system(&mut p, [10.0, -4.0], &cfg).expect("system path must still solve");
        assert!(
            sol.converged,
            "system path did not converge: {:?}",
            sol.reason
        );
        assert!(
            sol.n_cost_evals > 0,
            "system path spent no cost evaluations — the square-root flag is \
             fatal here rather than inert"
        );
        assert!((sol.x[0] - 1.0).abs() < 1e-8, "x0={}", sol.x[0]);
        assert!((sol.x[1] - 2.0).abs() < 1e-8, "x1={}", sol.x[1]);
    }

    #[test]
    fn test_zero_damping_diagonal() {
        // Zero Jacobian, no prior: μ-escalation is provably futile —
        // surfaced up-front on its own axis.
        let mut p = Tracked::new(|_x: &[f64; 2]| (vec![1.0], vec![[0.0, 0.0]]));
        let err = solve(&mut p, [0.0; 2], &config(), None).unwrap_err();
        assert!(
            matches!(err, LMError::ZeroDampingDiagonal { iteration: 0 }),
            "{err:?}"
        );
    }

    #[test]
    fn test_invalid_config_each_field() {
        let check = |cfg: LMConfig, want: fn(&ConfigDefect) -> bool| {
            let mut p_local = Tracked::new(|x: &[f64; 1]| (vec![x[0]], vec![[1.0]]));
            let err = solve(&mut p_local, [1.0], &cfg, None).unwrap_err();
            match err {
                LMError::InvalidConfig { defect } => assert!(want(&defect), "{defect:?}"),
                other => panic!("expected InvalidConfig, got {other:?}"),
            }
        };
        let mut c = config();
        c.max_iterations = 0;
        check(c, |d| matches!(d, ConfigDefect::MaxIterationsZero));
        let mut c = config();
        c.max_inner_trials = 0;
        check(c, |d| matches!(d, ConfigDefect::MaxInnerTrialsZero));
        let mut c = config();
        c.tau = 0.0;
        check(c, |d| matches!(d, ConfigDefect::TauNotPositive { .. }));
        let mut c = config();
        c.tau = f64::NAN;
        check(c, |d| matches!(d, ConfigDefect::TauNotPositive { .. }));
        let mut c = config();
        c.mu_max = 0.0;
        check(c, |d| matches!(d, ConfigDefect::MuMaxNotPositive { .. }));
        let mut c = config();
        c.min_relative_decrease = -1.0;
        check(c, |d| {
            matches!(d, ConfigDefect::MinRelativeDecreaseNegative { .. })
        });
        let mut c = config();
        c.gtol = -1.0;
        check(c, |d| matches!(d, ConfigDefect::GtolNegative { .. }));
        let mut c = config();
        c.qtol = -1.0;
        check(c, |d| matches!(d, ConfigDefect::QtolNegative { .. }));
        let mut c = config();
        c.xtol = f64::INFINITY;
        check(c, |d| matches!(d, ConfigDefect::XtolNegative { .. }));
        let mut c = config();
        c.ftol = -1.0;
        check(c, |d| matches!(d, ConfigDefect::FtolNegative { .. }));
        let mut c = config();
        c.max_consecutive_invalid = 0;
        check(c, |d| matches!(d, ConfigDefect::MaxConsecutiveInvalidZero));
    }

    #[test]
    fn test_invalid_prior() {
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0]], vec![[1.0]]));
        let prior = NLLSPrior {
            mean: [f64::NAN],
            covariance_inv: [[1.0]],
        };
        let err = solve(&mut p, [0.0], &config(), Some(&prior)).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidPrior {
                    defect: PriorDefect::NonFiniteMean { index: 0 }
                }
            ),
            "{err:?}"
        );

        let prior = NLLSPrior {
            mean: [0.0],
            covariance_inv: [[-1.0]],
        };
        let err = solve(&mut p, [0.0], &config(), Some(&prior)).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidPrior {
                    defect: PriorDefect::NegativeDiagonal { index: 0, .. }
                }
            ),
            "{err:?}"
        );
    }

    // ── Trial-failure recovery and persistence ──

    /// A propagation blowup at an aggressive trial step is a
    /// recoverable rejection, not a fit-killer.
    #[test]
    fn test_inf_trial_cost_recovers() {
        struct Walled {
            inner: Tracked<fn(&[f64; 1]) -> Resid1>,
        }
        impl CostProblem<1> for Walled {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                if x[0] > 50.0 {
                    return Ok(f64::INFINITY); // "propagation blew up"
                }
                self.inner.evaluate_cost(x)
            }
        }
        impl ResidualProblem<1> for Walled {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.inner.evaluate(x)
            }
        }
        // From x0 = −100 the log-overshoot GN step lands far past +50.
        fn log_overshoot(x: &[f64; 1]) -> (Vec<f64>, Vec<[f64; 1]>) {
            let d = x[0] - 5.0;
            let r = (1.0 + d * d).ln();
            let j = 2.0 * d / (1.0 + d * d);
            (vec![r], vec![[j]])
        }
        let mut p = Walled {
            inner: Tracked::new(log_overshoot as fn(&[f64; 1]) -> _),
        };
        let mut cfg = config();
        cfg.max_iterations = 300;
        // The far-side wall eats more than the default 5 consecutive
        // trials before damping shrinks the step under it.
        cfg.max_consecutive_invalid = 12;
        let sol = solve(&mut p, [-100.0], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 5.0).abs() < 1e-3, "x={}", sol.x[0]);
        assert!(sol.n_invalid_trials > 0, "wall should have been hit");
    }

    #[test]
    fn test_persistent_invalid_trials() {
        struct AlwaysFailsCost;
        impl CostProblem<1> for AlwaysFailsCost {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Err(TestError("always fails"))
            }
        }
        impl ResidualProblem<1> for AlwaysFailsCost {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 3.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 3.0) * (x[0] - 3.0),
                })
            }
        }
        let err = solve(&mut AlwaysFailsCost, [0.0], &config(), None).unwrap_err();
        match err {
            LMError::PersistentInvalidTrials {
                consecutive,
                last_source,
                best_x,
                best_cost,
                ..
            } => {
                assert_eq!(consecutive, 5);
                assert_eq!(last_source, Some(TestError("always fails")));
                assert_eq!(best_x, vec![0.0]);
                assert!((best_cost - 9.0).abs() < 1e-12);
            }
            other => panic!("expected PersistentInvalidTrials, got {other:?}"),
        }
    }

    // ── Termination reasons, one test each ──

    #[test]
    fn test_reason_gradient_tolerance() {
        // Linear: one exact step, then zero gradient at zero cost.
        let mut cfg = config();
        cfg.qtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0] - 3.0], vec![[1.0]]));
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();
        assert!(sol.converged);
        assert_eq!(sol.reason, TerminationReason::GradientTolerance);
        assert!((sol.x[0] - 3.0).abs() < 1e-10);
    }

    #[test]
    fn test_reason_step_tolerance_at_stationary_start() {
        // Start exactly at the minimum with the gradient test disabled:
        // the natural step is zero → StepTolerance with x unchanged.
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.ftol = 0.0;
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0] - 3.0], vec![[1.0]]));
        let sol = solve(&mut p, [3.0], &cfg, None).unwrap();
        assert!(sol.converged);
        assert_eq!(sol.reason, TerminationReason::StepTolerance);
        assert_eq!(sol.x[0], 3.0);
        assert!(sol.accepted_step_qnorm.is_none());
    }

    #[test]
    fn test_reason_max_iterations() {
        let mut cfg = config();
        cfg.max_iterations = 2;
        let mut p = Tracked::new(rosenbrock);
        let sol = solve(&mut p, [-1.2, 1.0], &cfg, None).unwrap();
        assert!(!sol.converged);
        assert_eq!(sol.reason, TerminationReason::MaxIterations);
        assert_eq!(sol.iterations, 2);
    }

    #[test]
    fn test_reason_damping_exhausted_on_adversarial_cost() {
        // The trial cost claims every step makes things worse while the
        // gradient says otherwise: every trial is rejected, μ escalates
        // to the cap, and the solver reports it — explicitly, finitely.
        struct Adversarial;
        impl CostProblem<1> for Adversarial {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(10.0)
            }
        }
        impl ResidualProblem<1> for Adversarial {
            fn evaluate(&mut self, _x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![1.0],
                    jacobian: vec![[1.0]],
                    cost: 1.0,
                })
            }
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.qtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        let sol = solve(&mut Adversarial, [5.0], &cfg, None).unwrap();
        assert!(!sol.converged);
        assert!(
            matches!(sol.reason, TerminationReason::DampingExhausted { .. }),
            "{:?}",
            sol.reason
        );
        assert_eq!(sol.x[0], 5.0, "iterate must not move on rejections");
    }

    #[test]
    fn test_reason_inner_trials_exhausted() {
        struct Adversarial;
        impl CostProblem<1> for Adversarial {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(10.0)
            }
        }
        impl ResidualProblem<1> for Adversarial {
            fn evaluate(&mut self, _x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![1.0],
                    jacobian: vec![[1.0]],
                    cost: 1.0,
                })
            }
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.qtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        cfg.max_inner_trials = 3;
        cfg.mu_max = 1e300;
        let sol = solve(&mut Adversarial, [5.0], &cfg, None).unwrap();
        assert!(!sol.converged);
        assert_eq!(
            sol.reason,
            TerminationReason::InnerTrialsExhausted { trials: 3 }
        );
    }

    #[test]
    fn test_reason_cost_tolerance() {
        // A genuine cost plateau: an irreducible constant residual
        // dominates Φ near the minimum, so the relative reductions —
        // including the UNDAMPED step's own predicted gain — fall
        // below ftol·Φ while the gradient and step tests are disabled.
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.qtol = 0.0;
        cfg.xtol = 0.0;
        cfg.max_iterations = 500;
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0] - 3.0, 10.0], vec![[1.0], [0.0]]));
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert_eq!(sol.reason, TerminationReason::CostTolerance);
        assert!((sol.x[0] - 3.0).abs() < 1e-2, "x={}", sol.x[0]);
    }

    // ── Clamps ──

    /// A clamp that shrinks every step below tolerance must NOT
    /// manufacture convergence (the old solver's documented false
    /// contract). With step reasons unreachable the solver runs its
    /// budget out while still making (slow) progress.
    #[test]
    fn test_clamped_step_cannot_declare_convergence() {
        struct Clamped {
            inner: Tracked<fn(&[f64; 1]) -> Resid1>,
        }
        impl CostProblem<1> for Clamped {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                self.inner.evaluate_cost(x)
            }
            fn constrain_step(&mut self, _x: &[f64; 1], delta: &mut [f64; 1]) {
                let cap = 0.01;
                if delta[0].abs() > cap {
                    delta[0] = if delta[0] > 0.0 { cap } else { -cap };
                }
            }
        }
        impl ResidualProblem<1> for Clamped {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.inner.evaluate(x)
            }
        }
        fn linear(x: &[f64; 1]) -> (Vec<f64>, Vec<[f64; 1]>) {
            (vec![x[0] - 1000.0], vec![[1.0]])
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.ftol = 0.0;
        cfg.xtol = 0.0;
        cfg.qtol = 1.0; // generous: the 0.01-clamped step is far below
        cfg.max_iterations = 50;
        let mut p = Clamped {
            inner: Tracked::new(linear as fn(&[f64; 1]) -> _),
        };
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();
        assert!(
            !sol.converged,
            "clamp manufactured convergence: {:?}",
            sol.reason
        );
        assert_eq!(sol.reason, TerminationReason::MaxIterations);
        // ... while still making real progress (one clamped accepted
        // step per outer iteration).
        assert!((sol.x[0] - 0.5).abs() < 1e-9, "x={}", sol.x[0]);
    }

    // ── Lifecycle and state discipline ──

    #[test]
    fn test_lifecycle_hooks_and_monotone_cost() {
        let mut p = Tracked::new(overshoot_rational);
        let sol = solve(&mut p, [6.5], &config(), None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 5.0).abs() < 1e-4);
        assert!(!p.rejected.is_empty(), "expected genuine rejections");
        assert!(!p.accepted.is_empty());
        // The final accepted hook saw the returned x.
        assert_eq!(p.accepted.last().unwrap()[0], sol.x[0]);
        // Full evaluations happened only at x0 + provisionally accepted
        // points; x0 itself is committed, so the counts match — and
        // every committed cost decreased monotonically.
        assert_eq!(p.full_evals.len(), p.accepted.len());
        assert_eq!(p.accepted[0][0], 6.5, "x0 is committed first");
        for w in p.full_evals.windows(2) {
            assert!(w[1] < w[0], "cost not monotone: {:?}", p.full_evals);
        }
        // No rejected trial point was ever committed.
        for r in &p.rejected {
            assert!(p.accepted.iter().all(|a| a != r));
        }
    }

    // ── Covariance ──

    #[test]
    fn test_covariance_known_linear() {
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0] - 3.0], vec![[1.0]]));
        let sol = solve(&mut p, [0.0], &config(), None).unwrap();
        let cov = sol.covariance.unwrap();
        assert!((cov[0][0] - 1.0).abs() < 1e-10);
    }

    /// Rank-deficient fit: converges in the observable subspace, keeps
    /// the solution, and reports the covariance failure explicitly —
    /// NEVER an all-zeros (σ = 0) matrix.
    #[test]
    fn test_covariance_singular_is_explicit_never_zeros() {
        let mut p = Tracked::new(|x: &[f64; 2]| (vec![x[0] - 3.0], vec![[1.0, 0.0]]));
        let sol = solve(&mut p, [0.0, 0.0], &config(), None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 3.0).abs() < 1e-8, "x0={}", sol.x[0]);
        assert_eq!(sol.x[1], 0.0);
        assert_eq!(
            sol.covariance.unwrap_err(),
            CovarianceFailure::SingularNormalMatrix
        );
    }

    // ── Mixed scales (equilibrated solve) ──

    #[test]
    fn test_extreme_unit_mixing() {
        // Parameter scales differ by ~10 orders of magnitude with a
        // coupling row — the equilibrated Cholesky must handle what a
        // raw-basis factorization cannot.
        let a = 3.0_f64;
        let b = 7e9_f64;
        let mut p = Tracked::new(move |x: &[f64; 2]| {
            (
                vec![
                    x[0] - a,
                    1e-10 * (x[1] - b),
                    1e-5 * (x[0] - a) + 1e-15 * (x[1] - b),
                ],
                vec![[1.0, 0.0], [0.0, 1e-10], [1e-5, 1e-15]],
            )
        });
        let mut cfg = config();
        cfg.max_iterations = 200;
        let sol = solve(&mut p, [0.0, 0.0], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - a).abs() < 1e-6, "x0={}", sol.x[0]);
        assert!((sol.x[1] - b).abs() / b < 1e-6, "x1={}", sol.x[1]);
        assert!(sol.covariance.is_ok());
    }

    // ── System path ──

    /// Quadratic objective Φ = (x−a)ᵀQ(x−a) + c via the
    /// normal-equations-level trait: A = Q, g = Q(a−x).
    struct QuadSystem {
        q: [[f64; 2]; 2],
        a: [f64; 2],
        c: f64,
    }
    impl QuadSystem {
        fn cost(&self, x: &[f64; 2]) -> f64 {
            let d = [x[0] - self.a[0], x[1] - self.a[1]];
            let mut phi = self.c;
            for i in 0..2 {
                for j in 0..2 {
                    phi += d[i] * self.q[i][j] * d[j];
                }
            }
            phi
        }
    }
    impl CostProblem<2> for QuadSystem {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 2]) -> Result<f64, TestError> {
            Ok(Self::cost(self, x))
        }
    }
    impl SystemProblem<2> for QuadSystem {
        fn evaluate_system(&mut self, x: &[f64; 2]) -> Result<SystemEvaluation<2>, TestError> {
            let mut rhs = [0.0; 2];
            for i in 0..2 {
                for j in 0..2 {
                    rhs[i] += self.q[i][j] * (self.a[j] - x[j]);
                }
            }
            Ok(SystemEvaluation {
                cost: Self::cost(self, x),
                normal: self.q,
                rhs,
            })
        }
    }

    #[test]
    fn test_solve_system_quadratic() {
        let mut p = QuadSystem {
            q: [[2.0, 0.0], [0.0, 8.0]],
            a: [1.0, 2.0],
            c: 3.0,
        };
        let sol = solve_system(&mut p, [10.0, -4.0], &config()).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 1.0).abs() < 1e-8, "x0={}", sol.x[0]);
        assert!((sol.x[1] - 2.0).abs() < 1e-8, "x1={}", sol.x[1]);
        assert!((sol.cost - 3.0).abs() < 1e-8);
        assert_eq!(sol.cost, sol.data_cost, "system path owns its objective");
        let cov = sol.covariance.unwrap();
        assert!((cov[0][0] - 0.5).abs() < 1e-10);
        assert!((cov[1][1] - 0.125).abs() < 1e-10);
    }

    #[test]
    fn test_invalid_system_negative_diagonal() {
        struct BadSystem;
        impl CostProblem<1> for BadSystem {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(1.0)
            }
        }
        impl SystemProblem<1> for BadSystem {
            fn evaluate_system(&mut self, _x: &[f64; 1]) -> Result<SystemEvaluation<1>, TestError> {
                Ok(SystemEvaluation {
                    cost: 1.0,
                    normal: [[-1.0]],
                    rhs: [0.0],
                })
            }
        }
        let err = solve_system(&mut BadSystem, [0.0], &config()).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidSystem {
                    iteration: 0,
                    defect: SystemDefect::NegativeDiagonal { index: 0, .. }
                }
            ),
            "{err:?}"
        );
    }

    #[test]
    fn test_invalid_system_nonfinite() {
        struct NanSystem;
        impl CostProblem<1> for NanSystem {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(1.0)
            }
        }
        impl SystemProblem<1> for NanSystem {
            fn evaluate_system(&mut self, _x: &[f64; 1]) -> Result<SystemEvaluation<1>, TestError> {
                Ok(SystemEvaluation {
                    cost: 1.0,
                    normal: [[f64::NAN]],
                    rhs: [0.0],
                })
            }
        }
        let err = solve_system(&mut NanSystem, [0.0], &config()).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidSystem {
                    iteration: 0,
                    defect: SystemDefect::NonFiniteNormal { row: 0, col: 0 }
                }
            ),
            "{err:?}"
        );
    }

    // ── Known-answer check on an exactly linear problem ──

    #[test]
    fn test_exact_solution_on_linear() {
        // Data lying exactly on \( y = 2x + 1 \): the least-squares
        // solution is \( p = [2, 1] \) with zero residual.
        let xs = [0.0, 1.0, 2.0, 3.0];
        let ys = [1.0, 3.0, 5.0, 7.0];

        let mut p = Tracked::new(move |p: &[f64; 2]| {
            let residuals: Vec<f64> = xs
                .iter()
                .zip(ys.iter())
                .map(|(&x, &y)| p[0] * x + p[1] - y)
                .collect();
            let jacobian: Vec<[f64; 2]> = xs.iter().map(|&x| [x, 1.0]).collect();
            (residuals, jacobian)
        });
        let mut cfg = config();
        cfg.xtol = 1e-14;
        let solution = solve(&mut p, [0.0; 2], &cfg, None).unwrap();

        assert!(solution.converged);
        assert!((solution.x[0] - 2.0).abs() < 1e-9);
        assert!((solution.x[1] - 1.0).abs() < 1e-9);
    }

    // ── Acceptance-test edge cases ──

    /// THE false-convergence regression that the 2020 CD3 capture-
    /// spanning fit exposed: on a stiff valley the rejections inflate
    /// μ by orders of magnitude, and a μ-crushed accepted step has a
    /// tiny quadratic form at an arbitrarily BAD iterate. Step-based
    /// convergence must therefore be judged on the UNDAMPED
    /// Gauss-Newton step — a μ-starved accepted step must never read
    /// as "converged".
    ///
    /// Adversarial cost: improvements exist only at micro-scale
    /// (|h| < 1e-7), so every natural-scale trial is rejected, μ
    /// inflates, and eventually micro-steps get accepted. The undamped
    /// GN step (and the gradient) remain large throughout — the solver
    /// must NOT declare convergence, no matter how generous qtol is.
    #[test]
    fn test_mu_starved_accepted_step_is_not_convergence() {
        // Positional cost: a microscopic basin at the origin inside a
        // wall — every macro trial is rejected (cost rises), only
        // |x| < 1e-7 micro-steps improve, while the residual/Jacobian
        // claim a huge gradient toward 1000 (so the undamped GN step
        // and the gradient stay enormous).
        fn stiff_cost(x: f64) -> f64 {
            if x.abs() < 1e-7 {
                1e6 - 1e-3 * (1e-7 - x.abs()) / 1e-7
            } else {
                1e6 + 1.0
            }
        }
        struct StiffValley;
        impl CostProblem<1> for StiffValley {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(stiff_cost(x[0]))
            }
        }
        impl ResidualProblem<1> for StiffValley {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 1000.0],
                    jacobian: vec![[1.0]],
                    cost: stiff_cost(x[0]),
                })
            }
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        cfg.qtol = 1.0; // generous: any micro-step is far below this
        cfg.max_iterations = 30;
        let sol = solve(&mut StiffValley, [0.0], &cfg, None).unwrap();
        // Whatever the budget outcome, a step-based "converged" at the
        // garbage iterate is forbidden.
        assert_ne!(sol.reason, TerminationReason::StepTolerance, "{sol:?}");
        assert_ne!(sol.reason, TerminationReason::CostTolerance, "{sol:?}");
        assert!(!sol.converged, "{:?}", sol.reason);
        assert!(sol.x[0].abs() < 1.0, "stayed near start: {}", sol.x[0]);
    }

    /// Monotonicity guard: when the trial cost claims an improvement
    /// but the full evaluation at the provisionally accepted point
    /// shows a HIGHER cost (trial/full objective inconsistency — the
    /// problem-contract violation the design doc warns about), the
    /// acceptance must be rolled back, nothing committed, and the
    /// inconsistency surfaced as persistent invalid trials.
    #[test]
    fn test_rollback_on_inconsistent_full_evaluation() {
        struct Inconsistent {
            full_calls: usize,
            committed: Vec<f64>,
        }
        impl CostProblem<1> for Inconsistent {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                // Honest, improving trial cost.
                Ok((x[0] - 3.0) * (x[0] - 3.0))
            }
            fn on_step_accepted(&mut self, x: &[f64; 1]) {
                self.committed.push(x[0]);
            }
        }
        impl ResidualProblem<1> for Inconsistent {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.full_calls += 1;
                if self.full_calls == 1 {
                    // Honest at x0.
                    return Ok(NLLSEvaluation {
                        residuals: vec![x[0] - 3.0],
                        jacobian: vec![[1.0]],
                        cost: (x[0] - 3.0) * (x[0] - 3.0),
                    });
                }
                // Full evaluation disagrees: claims things got WORSE.
                Ok(NLLSEvaluation {
                    residuals: vec![100.0],
                    jacobian: vec![[1.0]],
                    cost: 10000.0,
                })
            }
        }
        let mut p = Inconsistent {
            full_calls: 0,
            committed: Vec::new(),
        };
        let err = solve(&mut p, [0.0], &config(), None).unwrap_err();
        match err {
            LMError::PersistentInvalidTrials {
                best_x, best_cost, ..
            } => {
                assert_eq!(best_x, vec![0.0], "iterate must not move");
                assert!((best_cost - 9.0).abs() < 1e-12);
            }
            other => panic!("expected PersistentInvalidTrials, got {other:?}"),
        }
        assert_eq!(
            p.committed,
            vec![0.0],
            "only x0 may commit; rolled-back acceptances must not"
        );
    }

    /// pred ≤ 0 forces rejection even when the trial cost IMPROVES:
    /// accepting a step the local model calls bad would feed a
    /// negative ρ into Nielsen's update. A clamp that teleports the
    /// iterate to the mirror minimum (cost 0!) must still be rejected.
    #[test]
    fn test_pred_nonpositive_forces_rejection() {
        struct Teleport {
            inner: Tracked<fn(&[f64; 1]) -> Resid1>,
        }
        impl CostProblem<1> for Teleport {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                self.inner.evaluate_cost(x)
            }
            fn constrain_step(&mut self, x: &[f64; 1], delta: &mut [f64; 1]) {
                // Adversarial clamp: from the x = 3 region, jump
                // straight to the mirror minimum at −2.
                delta[0] = -2.0 - x[0];
            }
        }
        impl ResidualProblem<1> for Teleport {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.inner.evaluate(x)
            }
        }
        // r = x² − 4: minima at ±2. From x0 = 3 the model's descent
        // direction is toward +2; the clamped step to −2 has pred < 0.
        fn bimodal(x: &[f64; 1]) -> Resid1 {
            (vec![x[0] * x[0] - 4.0], vec![[2.0 * x[0]]])
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.ftol = 0.0;
        cfg.xtol = 0.0;
        cfg.qtol = 0.0;
        cfg.max_iterations = 3;
        let mut p = Teleport {
            inner: Tracked::new(bimodal as fn(&[f64; 1]) -> _),
        };
        let sol = solve(&mut p, [3.0], &cfg, None).unwrap();
        // Every clamped trial lands at −2 with cost 0 — a genuine
        // improvement — yet pred < 0 must reject every one of them.
        assert_eq!(sol.x[0], 3.0, "model-inconsistent step was accepted");
        assert!(!sol.converged);
        assert!(sol.n_rejected_trials > 0);
        assert!(p.inner.accepted.is_empty());
    }

    /// NaN trial costs follow the same invalid-trial path as INF.
    #[test]
    fn test_nan_trial_cost_is_invalid() {
        struct NanCost;
        impl CostProblem<1> for NanCost {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(f64::NAN)
            }
        }
        impl ResidualProblem<1> for NanCost {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 3.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 3.0) * (x[0] - 3.0),
                })
            }
        }
        let err = solve(&mut NanCost, [0.0], &config(), None).unwrap_err();
        match err {
            LMError::PersistentInvalidTrials {
                consecutive,
                last_source,
                ..
            } => {
                assert_eq!(consecutive, 5);
                assert_eq!(last_source, None, "NaN carries no domain error");
            }
            other => panic!("expected PersistentInvalidTrials, got {other:?}"),
        }
    }

    // ── Negative-cost validation (NaN-poisoned gradient regression) ──

    /// A negative cost would NaN-poison √Φ in the gradient test and
    /// read as instant false convergence — it must be a loud error on
    /// the residual path, including at cancellation scale.
    #[test]
    fn test_negative_cost_residual_path_is_error() {
        for bad in [-1.0, -1e-30] {
            struct NegCost(f64);
            impl CostProblem<1> for NegCost {
                type Error = TestError;
                fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                    Ok(self.0)
                }
            }
            impl ResidualProblem<1> for NegCost {
                fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                    Ok(NLLSEvaluation {
                        residuals: vec![x[0] - 1000.0],
                        jacobian: vec![[1.0]],
                        cost: self.0,
                    })
                }
            }
            let err = solve(&mut NegCost(bad), [0.0], &config(), None).unwrap_err();
            assert!(
                matches!(
                    err,
                    LMError::InvalidEvaluation {
                        iteration: 0,
                        defect: EvaluationDefect::NegativeCost { .. }
                    }
                ),
                "cost={bad}: {err:?}"
            );
        }
    }

    /// Same regression on the system path — the realistic trigger is
    /// a Schur profiled cost rounding negative through cancellation.
    #[test]
    fn test_negative_cost_system_path_is_error() {
        for bad in [-5.0, -1e-30] {
            struct NegSystem(f64);
            impl CostProblem<1> for NegSystem {
                type Error = TestError;
                fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                    Ok(self.0)
                }
            }
            impl SystemProblem<1> for NegSystem {
                fn evaluate_system(
                    &mut self,
                    _x: &[f64; 1],
                ) -> Result<SystemEvaluation<1>, TestError> {
                    Ok(SystemEvaluation {
                        cost: self.0,
                        normal: [[1.0]],
                        rhs: [1000.0],
                    })
                }
            }
            let err = solve_system(&mut NegSystem(bad), [0.0], &config()).unwrap_err();
            assert!(
                matches!(
                    err,
                    LMError::InvalidSystem {
                        iteration: 0,
                        defect: SystemDefect::NegativeCost { .. }
                    }
                ),
                "cost={bad}: {err:?}"
            );
        }
    }

    // ── Per-variant coverage for the remaining defect axes ──

    #[test]
    fn test_nonfinite_residual_and_jacobian_entries_are_errors() {
        // The cost must be finite so the per-entry checks (not the
        // cost check) are what fire.
        struct NanResidual;
        impl CostProblem<1> for NanResidual {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(1.0)
            }
        }
        impl ResidualProblem<1> for NanResidual {
            fn evaluate(&mut self, _x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![f64::NAN],
                    jacobian: vec![[0.0]],
                    cost: 1.0,
                })
            }
        }
        let err = solve(&mut NanResidual, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidEvaluation {
                    iteration: 0,
                    defect: EvaluationDefect::Residual { index: 0 }
                }
            ),
            "{err:?}"
        );

        let mut p = Tracked::new(|_x: &[f64; 1]| (vec![1.0], vec![[f64::INFINITY]]));
        let err = solve(&mut p, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidEvaluation {
                    iteration: 0,
                    defect: EvaluationDefect::JacobianEntry { row: 0, col: 0 }
                }
            ),
            "{err:?}"
        );
    }

    #[test]
    fn test_system_nonfinite_cost_and_rhs_are_errors() {
        struct BadSys {
            cost: f64,
            rhs: f64,
        }
        impl CostProblem<1> for BadSys {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(1.0)
            }
        }
        impl SystemProblem<1> for BadSys {
            fn evaluate_system(&mut self, _x: &[f64; 1]) -> Result<SystemEvaluation<1>, TestError> {
                Ok(SystemEvaluation {
                    cost: self.cost,
                    normal: [[1.0]],
                    rhs: [self.rhs],
                })
            }
        }
        let err = solve_system(
            &mut BadSys {
                cost: f64::NAN,
                rhs: 0.0,
            },
            [0.0],
            &config(),
        )
        .unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidSystem {
                    iteration: 0,
                    defect: SystemDefect::NonFiniteCost
                }
            ),
            "{err:?}"
        );

        let err = solve_system(
            &mut BadSys {
                cost: 1.0,
                rhs: f64::NAN,
            },
            [0.0],
            &config(),
        )
        .unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidSystem {
                    iteration: 0,
                    defect: SystemDefect::NonFiniteRhs { index: 0 }
                }
            ),
            "{err:?}"
        );
    }

    #[test]
    fn test_invalid_prior_nonfinite_covariance_inv() {
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0]], vec![[1.0]]));
        let prior = NLLSPrior {
            mean: [0.0],
            covariance_inv: [[f64::NAN]],
        };
        let err = solve(&mut p, [0.0], &config(), Some(&prior)).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidPrior {
                    defect: PriorDefect::NonFiniteCovarianceInv { row: 0, col: 0 }
                }
            ),
            "{err:?}"
        );
    }

    #[test]
    fn test_invalid_config_tau_exceeds_mu_max() {
        let mut cfg = config();
        cfg.tau = 1.0;
        cfg.mu_max = 0.5;
        let mut p = Tracked::new(|x: &[f64; 1]| (vec![x[0]], vec![[1.0]]));
        let err = solve(&mut p, [1.0], &cfg, None).unwrap_err();
        assert!(
            matches!(
                err,
                LMError::InvalidConfig {
                    defect: ConfigDefect::TauExceedsMuMax { .. }
                }
            ),
            "{err:?}"
        );
    }

    /// An asymmetric (at rounding level) prior precision matrix — the
    /// realistic LU-inversion artifact — is symmetrized on receipt,
    /// not rejected, and the solve proceeds to the exact MAP point of
    /// the symmetrized objective.
    #[test]
    fn test_asymmetric_prior_is_symmetrized_not_rejected() {
        let mut p = Tracked::new(|x: &[f64; 2]| {
            (vec![x[0] - 10.0, x[1] - 10.0], vec![[1.0, 0.0], [0.0, 1.0]])
        });
        // Off-diagonal asymmetric by 1e-16-scale, as an LU-derived
        // inverse would be.
        let prior = NLLSPrior {
            mean: [0.0, 0.0],
            covariance_inv: [[1.0, 0.1 + 1e-16], [0.1 - 1e-16, 1.0]],
        };
        let sol = solve(&mut p, [0.0, 0.0], &config(), Some(&prior)).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        // MAP of (x−10)² + xᵀPx per component with symmetrized
        // P = [[1, 0.1], [0.1, 1]]: solve (I + P) x = (10, 10).
        // (2, 0.1; 0.1, 2) x = (10,10) → x = 10/2.1 each.
        let expected = 10.0 / 2.1;
        assert!((sol.x[0] - expected).abs() < 1e-6, "x0={}", sol.x[0]);
        assert!((sol.x[1] - expected).abs() < 1e-6, "x1={}", sol.x[1]);
    }

    // ── Rollback axes beyond cost inconsistency ──

    /// Domain error from the full evaluation at a provisionally
    /// accepted point: rolled back, surfaced as persistent invalid
    /// trials with the domain error retained.
    #[test]
    fn test_rollback_on_domain_error_at_accepted_point() {
        struct FailsAfterFirst {
            full_calls: usize,
        }
        impl CostProblem<1> for FailsAfterFirst {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok((x[0] - 3.0) * (x[0] - 3.0))
            }
        }
        impl ResidualProblem<1> for FailsAfterFirst {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.full_calls += 1;
                if self.full_calls > 1 {
                    return Err(TestError("propagation died at accepted point"));
                }
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 3.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 3.0) * (x[0] - 3.0),
                })
            }
        }
        let err = solve(
            &mut FailsAfterFirst { full_calls: 0 },
            [0.0],
            &config(),
            None,
        )
        .unwrap_err();
        match err {
            LMError::PersistentInvalidTrials {
                last_source,
                best_x,
                ..
            } => {
                assert_eq!(
                    last_source,
                    Some(TestError("propagation died at accepted point"))
                );
                assert_eq!(best_x, vec![0.0]);
            }
            other => panic!("expected PersistentInvalidTrials, got {other:?}"),
        }
    }

    /// A mid-solve CONTRACT violation (dimension mismatch at an
    /// accepted point) is fatal on its own axis with the right
    /// iteration — not recycled into a rollback.
    #[test]
    fn test_hard_defect_mid_solve_is_fatal() {
        struct BreaksContract {
            full_calls: usize,
        }
        impl CostProblem<1> for BreaksContract {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok((x[0] - 3.0) * (x[0] - 3.0))
            }
        }
        impl ResidualProblem<1> for BreaksContract {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.full_calls += 1;
                if self.full_calls > 1 {
                    return Ok(NLLSEvaluation {
                        residuals: vec![1.0, 2.0],
                        jacobian: vec![[1.0]],
                        cost: 5.0,
                    });
                }
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 3.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 3.0) * (x[0] - 3.0),
                })
            }
        }
        let err = solve(
            &mut BreaksContract { full_calls: 0 },
            [0.0],
            &config(),
            None,
        )
        .unwrap_err();
        assert!(
            matches!(
                err,
                LMError::DimensionMismatch {
                    iteration: 1,
                    residuals: 2,
                    jacobian: 1
                }
            ),
            "{err:?}"
        );
    }

    /// An overflowing predicted reduction (INF/NaN from an adversarial
    /// clamp) is a forced rejection — never a panic, never an accept.
    #[test]
    fn test_pred_overflow_forces_rejection() {
        struct HugeClamp {
            inner: Tracked<fn(&[f64; 1]) -> Resid1>,
        }
        impl CostProblem<1> for HugeClamp {
            type Error = TestError;
            fn evaluate_cost(&mut self, _x: &[f64; 1]) -> Result<f64, TestError> {
                // Adversarial liar: claims a perfect fit at the
                // teleported point, so ONLY the pred > 0 guard stands
                // between the solver and a model-inconsistent accept.
                Ok(0.0)
            }
            fn constrain_step(&mut self, _x: &[f64; 1], delta: &mut [f64; 1]) {
                delta[0] = 1e200;
            }
        }
        impl ResidualProblem<1> for HugeClamp {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.inner.evaluate(x)
            }
        }
        fn scaled(x: &[f64; 1]) -> Resid1 {
            (vec![1e150 * (x[0] - 1.0)], vec![[1e150]])
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.qtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        cfg.max_iterations = 3;
        let mut p = HugeClamp {
            inner: Tracked::new(scaled as fn(&[f64; 1]) -> _),
        };
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();
        assert_eq!(sol.x[0], 0.0, "overflowed-pred step must not be accepted");
        assert!(!sol.converged);
    }

    /// Problem-side committed state is BIT-identical across rejected
    /// trials: stage in evaluate_cost, commit only in
    /// on_step_accepted, and verify the committed bits never moved on
    /// an always-rejecting run.
    #[test]
    fn test_committed_state_bit_identical_across_rejections() {
        struct Staged {
            committed: u64,
            pending: u64,
        }
        impl CostProblem<1> for Staged {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                self.pending = x[0].to_bits();
                Ok(10.0) // adversarial: every trial is rejected
            }
            fn on_step_accepted(&mut self, _x: &[f64; 1]) {
                self.committed = self.pending;
            }
        }
        impl ResidualProblem<1> for Staged {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                self.pending = x[0].to_bits();
                Ok(NLLSEvaluation {
                    residuals: vec![1.0],
                    jacobian: vec![[1.0]],
                    cost: 1.0,
                })
            }
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.qtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        let mut p = Staged {
            committed: u64::MAX,
            pending: u64::MAX,
        };
        let x0 = 5.0_f64;
        let sol = solve(&mut p, [x0], &cfg, None).unwrap();
        assert!(!sol.converged);
        assert_eq!(
            p.committed,
            x0.to_bits(),
            "committed state must still be the x0 commit, bit-for-bit"
        );
    }

    // ── Geodesic acceleration ──

    /// Rosenbrock with the EXACT directional second derivative
    /// supplied by the problem (what a one-parameter Jet2 evaluation
    /// provides in production): r₁ = 10(x₂ − x₁²) has
    /// ∂²r₁/∂x₁² = −20, so r''_vv = [−20 v₁², 0].
    struct RosenbrockGeo {
        hook_calls: std::cell::Cell<usize>,
    }
    impl CostProblem<2> for RosenbrockGeo {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 2]) -> Result<f64, TestError> {
            let r1 = 10.0 * (x[1] - x[0] * x[0]);
            let r2 = 1.0 - x[0];
            Ok(r1 * r1 + r2 * r2)
        }
        fn second_directional_derivative(
            &mut self,
            _x: &[f64; 2],
            v: &[f64; 2],
        ) -> Option<Vec<f64>> {
            self.hook_calls.set(self.hook_calls.get() + 1);
            Some(vec![-20.0 * v[0] * v[0], 0.0])
        }
    }
    impl ResidualProblem<2> for RosenbrockGeo {
        fn evaluate(&mut self, x: &[f64; 2]) -> Result<NLLSEvaluation<2>, TestError> {
            let r1 = 10.0 * (x[1] - x[0] * x[0]);
            let r2 = 1.0 - x[0];
            Ok(NLLSEvaluation {
                residuals: vec![r1, r2],
                jacobian: vec![[-20.0 * x[0], 10.0], [-1.0, 0.0]],
                cost: r1 * r1 + r2 * r2,
            })
        }
    }

    /// Acceleration must reach the same minimum in fewer cost
    /// evaluations on the canonical curved-valley problem.
    #[test]
    fn test_geodesic_accelerates_rosenbrock() {
        let run = |geodesic: bool| {
            let mut cfg = config();
            cfg.geodesic_acceleration = geodesic;
            cfg.max_iterations = 500;
            let mut p = RosenbrockGeo {
                hook_calls: std::cell::Cell::new(0),
            };
            let sol = solve(&mut p, [-1.2, 1.0], &cfg, None).unwrap();
            (sol, p.hook_calls.get())
        };
        let (off, off_calls) = run(false);
        let (on, on_calls) = run(true);
        assert_eq!(off_calls, 0, "hook must not be called with the flag off");
        assert!(on_calls > 0, "hook must be exercised with the flag on");
        assert!(off.converged && on.converged);
        assert!((on.x[0] - 1.0).abs() < 1e-6 && (on.x[1] - 1.0).abs() < 1e-6);
        assert!((off.x[0] - 1.0).abs() < 1e-6 && (off.x[1] - 1.0).abs() < 1e-6);
        assert!(on.n_accelerated_trials > 0);
        assert!(
            on.n_cost_evals < off.n_cost_evals,
            "acceleration must reduce cost evaluations: on={} off={}",
            on.n_cost_evals,
            off.n_cost_evals
        );
    }

    /// Geodesic acceleration must survive the SQUARE-ROOT path.
    ///
    /// The acceleration solve reuses the Cholesky factor of the damped
    /// system; the square-root path never forms one by squaring, so it
    /// hands back R^T from the QR instead (R^T R = the scaled damped
    /// normal matrix). This pins that the hand-off works end to end:
    /// acceleration fires, is applied, and the solve converges to the
    /// same answer — the regression here would be silent (acceleration
    /// quietly declining on every trial), so the assertion is on the
    /// applied-trial COUNT, not on convergence alone.
    #[test]
    fn test_geodesic_acceleration_fires_on_the_square_root_path() {
        let mut cfg = config();
        cfg.geodesic_acceleration = true;
        cfg.square_root_solve = true;
        cfg.max_iterations = 500;
        let mut p = RosenbrockGeo {
            hook_calls: std::cell::Cell::new(0),
        };
        let sol = solve(&mut p, [-1.2, 1.0], &cfg, None).unwrap();
        assert!(sol.converged, "sqrt-path Rosenbrock must converge");
        assert!(
            (sol.x[0] - 1.0).abs() < 1e-6 && (sol.x[1] - 1.0).abs() < 1e-6,
            "sqrt-path answer moved: {:?}",
            sol.x
        );
        assert!(
            p.hook_calls.get() > 0,
            "second-derivative hook never exercised on the sqrt path"
        );
        assert!(
            sol.n_accelerated_trials > 0,
            "acceleration computed but never APPLIED on the sqrt path — the \
             R^T Cholesky hand-off has regressed to silently declining"
        );
    }

    /// The factor the square-root path hands to the acceleration solve
    /// must be the LOWER-triangular \\(R^\top\\), and solving through it
    /// must reproduce the normal-equations solve exactly.
    ///
    /// # Why this is asserted numerically and not through a fit
    ///
    /// [`solve_with_factor`] reads its argument as a lower-triangular
    /// Cholesky factor: below-diagonal entries in the forward substitution,
    /// above-diagonal in the back substitution. Hand it \\(R\\) instead of
    /// \\(R^\top\\) and every below-diagonal read lands on a structural
    /// zero, so the triangular solve silently degenerates into a DIAGONAL
    /// one — a finite, plausible, wrong acceleration. Nothing throws.
    ///
    /// The end-to-end acceleration gate does not catch that: Levenberg-
    /// Marquardt is robust to a poor acceleration, because the velocity
    /// step dominates and the gain-ratio test rejects the trials the
    /// correction spoils. The solve still converges and
    /// `n_accelerated_trials` is still positive. Only a direct comparison
    /// against the other path's solve can see it, so that is what this
    /// does — one damped system, both routes to \\(a\\), asserted equal.
    #[test]
    fn test_square_root_factor_reproduces_the_cholesky_solve() {
        const N: usize = 4;
        // A system with genuinely coupled columns, so a degenerate
        // "diagonal only" solve cannot coincidentally agree.
        let rows: Vec<[f64; N]> = vec![
            [1.0, 0.5, -0.25, 0.125],
            [0.0, 2.0, 0.75, -0.5],
            [-1.5, 0.25, 3.0, 1.0],
            [0.5, -1.0, 0.5, 2.5],
            [2.0, 1.0, -1.0, 0.25],
        ];
        let rhs = [1.0_f64, -2.0, 0.5, 3.0, -1.5];

        let mut normal = [[0.0_f64; N]; N];
        let mut grad = [0.0_f64; N];
        for (row, b_i) in rows.iter().zip(&rhs) {
            for j in 0..N {
                for k in 0..N {
                    normal[j][k] += row[j] * row[k];
                }
                grad[j] += row[j] * b_i;
            }
        }
        let normal = mat_symmetrize(&normal);
        let mut d = [0.0_f64; N];
        let d_max = update_scaling(&mut d, &normal);

        // An arbitrary right-hand side standing in for the acceleration
        // solve's -Jᵀ r''_vv, which has nothing to do with `grad`.
        let rhs_a = [0.3_f64, -1.7, 2.2, 0.9];

        for mu in [0.0_f64, 1e-3, 1.0, 1e6] {
            let mut qr = crate::linalg::qr::QrAccumulator::<N>::new();
            for (row, b_i) in rows.iter().zip(&rhs) {
                let mut scaled = [0.0_f64; N];
                for j in 0..N {
                    scaled[j] = row[j] / effective_scale(d[j], d_max);
                }
                qr.push_row(&scaled, *b_i);
            }
            qr.push_damping(mu);

            // The hand-off exactly as `solve_core` performs it.
            let r = qr.r();
            let mut l = [[0.0_f64; N]; N];
            for (i, row) in l.iter_mut().enumerate() {
                for (j, slot) in row.iter_mut().enumerate() {
                    *slot = r[j][i];
                }
            }
            let via_qr = solve_with_factor(&l, &rhs_a, &d, d_max).expect("R^T is a valid factor");
            let via_cholesky =
                solve_damped(&normal, &rhs_a, &d, d_max, mu).expect("damped system is definite");

            for i in 0..N {
                let scale = via_cholesky[i].abs().max(1e-12);
                assert!(
                    (via_qr[i] - via_cholesky[i]).abs() <= 1e-9 * scale,
                    "mu {mu:.0e} component {i}: R^T route gave {} but the Cholesky route \
                     gave {} — the factor handed to the acceleration solve is not R^T",
                    via_qr[i],
                    via_cholesky[i],
                );
            }
        }
    }

    /// An untrustworthy expansion (huge curvature) trips the avmax
    /// guard: the trial is rejected through μ escalation with no cost
    /// spent, no acceleration is ever applied, and the solve still
    /// finishes honestly.
    #[test]
    fn test_geodesic_avmax_guard_rejects_huge_curvature() {
        struct HugeCurvature;
        impl CostProblem<1> for HugeCurvature {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok((x[0] - 3.0) * (x[0] - 3.0))
            }
            fn second_directional_derivative(
                &mut self,
                _x: &[f64; 1],
                _v: &[f64; 1],
            ) -> Option<Vec<f64>> {
                Some(vec![1e12])
            }
        }
        impl ResidualProblem<1> for HugeCurvature {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 3.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 3.0) * (x[0] - 3.0),
                })
            }
        }
        let mut cfg = config();
        cfg.geodesic_acceleration = true;
        let sol = solve(&mut HugeCurvature, [0.0], &cfg, None).unwrap();
        assert_eq!(sol.n_accelerated_trials, 0);
        assert!(sol.n_rejected_trials > 0, "avmax violations must reject");
        // The guard never blocks honest termination.
        assert!(!sol.converged || (sol.x[0] - 3.0).abs() < 1e-6, "{sol:?}");
    }

    /// A hook returning the wrong number of rows is a contract
    /// violation — loud, on its own axis.
    #[test]
    fn test_geodesic_wrong_length_is_error() {
        struct WrongLen;
        impl CostProblem<1> for WrongLen {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok((x[0] - 3.0) * (x[0] - 3.0))
            }
            fn second_directional_derivative(
                &mut self,
                _x: &[f64; 1],
                _v: &[f64; 1],
            ) -> Option<Vec<f64>> {
                Some(vec![0.0, 0.0, 0.0])
            }
        }
        impl ResidualProblem<1> for WrongLen {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 3.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 3.0) * (x[0] - 3.0),
                })
            }
        }
        let mut cfg = config();
        cfg.geodesic_acceleration = true;
        let err = solve(&mut WrongLen, [0.0], &cfg, None).unwrap_err();
        assert!(matches!(err, LMError::DimensionMismatch { .. }), "{err:?}");
    }

    /// The normal-equations path carries no Jacobian rows, so the
    /// hook must never fire there even with the flag on.
    #[test]
    fn test_geodesic_inert_on_system_path() {
        struct CountingQuad {
            inner: QuadSystem,
            hook_calls: usize,
        }
        impl CostProblem<2> for CountingQuad {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 2]) -> Result<f64, TestError> {
                self.inner.evaluate_cost(x)
            }
            fn second_directional_derivative(
                &mut self,
                _x: &[f64; 2],
                _v: &[f64; 2],
            ) -> Option<Vec<f64>> {
                self.hook_calls += 1;
                Some(vec![0.0])
            }
        }
        impl SystemProblem<2> for CountingQuad {
            fn evaluate_system(&mut self, x: &[f64; 2]) -> Result<SystemEvaluation<2>, TestError> {
                self.inner.evaluate_system(x)
            }
        }
        let mut cfg = config();
        cfg.geodesic_acceleration = true;
        let mut p = CountingQuad {
            inner: QuadSystem {
                q: [[2.0, 0.0], [0.0, 8.0]],
                a: [1.0, 2.0],
                c: 3.0,
            },
            hook_calls: 0,
        };
        let sol = solve_system(&mut p, [10.0, -4.0], &cfg).unwrap();
        assert!(sol.converged);
        assert_eq!(p.hook_calls, 0, "system path must never call the hook");
        assert_eq!(sol.n_accelerated_trials, 0);
    }

    // ── Determinism ──

    /// Golden-trace bit-determinism: the driver is fixed-order scalar
    /// f64 arithmetic (no libm), so on any IEEE 754 platform this
    /// rational-arithmetic problem must reproduce these exact bits.
    /// A failure here means the driver's arithmetic or its evaluation
    /// ORDER changed — treat as a breaking change, not a flake.
    #[test]
    fn test_golden_trace_bit_determinism() {
        let run = || {
            let mut p = Tracked::new(overshoot_rational);
            solve(&mut p, [6.5], &config(), None).unwrap()
        };
        let a = run();
        let b = run();
        assert_eq!(a.x[0].to_bits(), b.x[0].to_bits());
        assert_eq!(a.cost.to_bits(), b.cost.to_bits());
        assert_eq!(a.mu_final.to_bits(), b.mu_final.to_bits());
        assert_eq!(a.iterations, b.iterations);
        assert_eq!(a.n_cost_evals, b.n_cost_evals);
        assert_eq!(a.n_rejected_trials, b.n_rejected_trials);

        // Cross-platform pins (generated on macOS arm64; must hold on
        // every IEEE 754 platform — see module docs on determinism).
        assert_eq!(a.x[0].to_bits(), GOLDEN_X_BITS, "x = {:?}", a.x[0]);
        assert_eq!(a.cost.to_bits(), GOLDEN_COST_BITS, "cost = {:?}", a.cost);
        assert_eq!(
            a.mu_final.to_bits(),
            GOLDEN_MU_FINAL_BITS,
            "mu = {:?}",
            a.mu_final
        );
        assert_eq!(a.iterations, GOLDEN_ITERATIONS);
        assert_eq!(a.n_cost_evals, GOLDEN_COST_EVALS);
        assert_eq!(a.n_rejected_trials, GOLDEN_REJECTED);

        // The full committed-iterate SEQUENCE is part of the trace —
        // not just the endpoint — so an evaluation-order change
        // anywhere in the driver shows up here.
        let mut p = Tracked::new(overshoot_rational);
        let _ = solve(&mut p, [6.5], &config(), None).unwrap();
        let accepted_bits: Vec<u64> = p.accepted.iter().map(|v| v[0].to_bits()).collect();
        let rejected_bits: Vec<u64> = p.rejected.iter().map(|v| v[0].to_bits()).collect();
        assert_eq!(
            accepted_bits, GOLDEN_ACCEPTED_BITS,
            "accepted iterate sequence changed"
        );
        assert_eq!(
            rejected_bits, GOLDEN_REJECTED_BITS,
            "rejected trial sequence changed"
        );
    }

    // ── Model refresh (inexact-model problems) ──

    /// A problem whose objective is an **anchored surrogate**: the
    /// residual is measured against a reference the problem is allowed
    /// to move — but only inside `refresh_model`, never per call.
    /// Every callback is logged so the driver's refresh sequence can
    /// be asserted exactly.
    struct AnchoredSurrogate {
        /// The anchor the current objective is measured against.
        anchor: f64,
        /// Anchor to adopt at the next refresh (one-shot).
        next_anchor: Option<f64>,
        /// `Some` to make the re-assembly after the refresh fail.
        fail_after_refresh: bool,
        log: Vec<(&'static str, f64)>,
    }

    impl AnchoredSurrogate {
        fn new(anchor: f64) -> Self {
            Self {
                anchor,
                next_anchor: None,
                fail_after_refresh: false,
                log: Vec::new(),
            }
        }
        fn residual(&self, x: f64) -> f64 {
            x - self.anchor
        }
        fn steps(&self, kind: &str) -> Vec<f64> {
            self.log
                .iter()
                .filter(|(k, _)| *k == kind)
                .map(|(_, v)| *v)
                .collect()
        }
    }

    impl CostProblem<1> for AnchoredSurrogate {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
            self.log.push(("cost", x[0]));
            let r = self.residual(x[0]);
            Ok(r * r)
        }
        fn on_step_accepted(&mut self, x: &[f64; 1]) {
            self.log.push(("accepted", x[0]));
        }
        fn on_step_rejected(&mut self, x: &[f64; 1]) {
            self.log.push(("rejected", x[0]));
        }
        fn refresh_model(&mut self, x: &[f64; 1]) -> bool {
            self.log.push(("refresh", x[0]));
            match self.next_anchor.take() {
                Some(a) => {
                    self.anchor = a;
                    true
                }
                None => false,
            }
        }
    }

    impl ResidualProblem<1> for AnchoredSurrogate {
        fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
            self.log.push(("evaluate", x[0]));
            if self.fail_after_refresh && self.next_anchor.is_none() && self.anchor != 1.0 {
                return Err(TestError("refreshed model cannot be evaluated"));
            }
            let r = self.residual(x[0]);
            Ok(NLLSEvaluation {
                residuals: vec![r],
                jacobian: vec![[1.0]],
                cost: r * r,
            })
        }
    }

    /// A `refresh_model` declaration must re-assemble at the CURRENT
    /// accepted point and re-commit it BEFORE any trial of that
    /// iteration — otherwise the driver would damp a system built from
    /// the old model while comparing trials evaluated under the new
    /// one, which is exactly the cross-model gain ratio the hook
    /// exists to prevent.
    #[test]
    fn test_refresh_model_reassembles_and_recommits_before_any_trial() {
        let mut p = AnchoredSurrogate::new(1.0);
        p.next_anchor = Some(3.0);
        let sol = solve(&mut p, [0.0], &config(), None).unwrap();

        // Refresh happened once, at the accepted point x = 0.
        assert_eq!(p.steps("refresh")[0], 0.0);
        assert_eq!(sol.n_model_refreshes, 1, "one declaration, one assembly");
        // ... and was immediately followed by a full evaluation and a
        // commit at that same point, before the first trial cost.
        let order: Vec<&'static str> = p.log.iter().map(|(k, _)| *k).collect();
        let first_refresh = order.iter().position(|k| *k == "refresh").unwrap();
        assert_eq!(
            &order[first_refresh..first_refresh + 3],
            &["refresh", "evaluate", "accepted"],
            "log: {:?}",
            p.log
        );
        assert_eq!(
            p.log[first_refresh + 1].1,
            0.0,
            "re-assembled at x, not x+h"
        );
        assert_eq!(p.log[first_refresh + 2].1, 0.0, "re-committed at x");

        // The solve minimized the REFRESHED objective, not the one it
        // started with.
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 3.0).abs() < 1e-6, "x = {}", sol.x[0]);
        assert!(sol.cost < 1e-12, "cost = {}", sol.cost);
    }

    /// The default hook is inert: a problem that never refreshes must
    /// spend no extra evaluation and must trace exactly as before.
    #[test]
    fn test_refresh_model_false_spends_no_evaluation() {
        let mut p = AnchoredSurrogate::new(3.0);
        let sol = solve(&mut p, [0.0], &config(), None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        // One evaluation at x0 plus one per accepted trial — no
        // re-assembly anywhere.
        let evaluates = p.steps("evaluate");
        let costs = p.steps("cost");
        assert_eq!(
            evaluates.len(),
            costs.len() + 1,
            "an extra full evaluation appeared without a refresh: {:?}",
            p.log
        );
        assert!(!p.steps("refresh").is_empty(), "hook was never called");
        assert_eq!(sol.n_model_refreshes, 0, "nothing declared a refresh");
    }

    /// The evaluation budget documented on
    /// [`LMConfig::max_iterations`] must match what the solver spends:
    /// a refresh costs one FULL assembly (residuals + Jacobian) beyond
    /// the trial-driven ones, and `LMSolution::n_model_refreshes` is
    /// the only field that reports it.
    ///
    /// Measured against a NO-OP refresh — the hook returns `true` while
    /// re-anchoring the surrogate at the value it already had. The
    /// objective is bit-identical before and after, and the refresh
    /// fires at iteration 1 where `last_accepted` is already `None`, so
    /// the post-refresh convergence battery re-runs the test that just
    /// returned `None` and the trajectory is unchanged. Anything the
    /// two solves differ by is therefore the refresh's own cost, and
    /// nothing else — which pins the count exactly at one extra
    /// assembly per declaration.
    #[test]
    fn test_refresh_model_assembly_is_counted_and_costs_one_evaluation() {
        // Baseline: same problem, hook never declares.
        let mut plain = AnchoredSurrogate::new(3.0);
        let sol_plain = solve(&mut plain, [0.0], &config(), None).unwrap();

        // Same problem, one no-op refresh declared at iteration 1.
        let mut refreshed = AnchoredSurrogate::new(3.0);
        refreshed.next_anchor = Some(3.0);
        let sol_refresh = solve(&mut refreshed, [0.0], &config(), None).unwrap();

        // The refresh changed nothing about the objective, so every
        // other counter and the whole trajectory must be identical.
        assert_eq!(sol_refresh.x[0].to_bits(), sol_plain.x[0].to_bits());
        assert_eq!(sol_refresh.cost.to_bits(), sol_plain.cost.to_bits());
        assert_eq!(sol_refresh.iterations, sol_plain.iterations);
        assert_eq!(sol_refresh.n_cost_evals, sol_plain.n_cost_evals);
        assert_eq!(sol_refresh.n_rejected_trials, sol_plain.n_rejected_trials);
        assert_eq!(sol_refresh.n_invalid_trials, sol_plain.n_invalid_trials);
        assert_eq!(sol_refresh.reason, sol_plain.reason);

        // The declaration is counted...
        assert_eq!(sol_plain.n_model_refreshes, 0);
        assert_eq!(sol_refresh.n_model_refreshes, 1);

        // ... and costs EXACTLY that many extra full assemblies.
        let plain_assemblies = plain.steps("evaluate").len();
        let refresh_assemblies = refreshed.steps("evaluate").len();
        assert_eq!(
            refresh_assemblies,
            plain_assemblies + sol_refresh.n_model_refreshes,
            "refresh assemblies {refresh_assemblies} vs baseline {plain_assemblies}, \
             reported refreshes {}; log: {:?}",
            sol_refresh.n_model_refreshes,
            refreshed.log,
        );

        // And the full documented budget line closes: assemblies =
        // 1 initial + provisional acceptances + refreshes. This
        // problem is exactly linear, so no acceptance is ever rolled
        // back and each accepted trial assembles once; the `accepted`
        // callbacks are the initial commit, the refresh re-commit, and
        // one per committed step.
        let commits = refreshed.steps("accepted").len() - 1 - sol_refresh.n_model_refreshes;
        assert_eq!(
            refresh_assemblies,
            1 + commits + sol_refresh.n_model_refreshes,
            "log: {:?}",
            refreshed.log
        );
        // Bounded by the iteration count, as documented.
        assert!(sol_refresh.n_model_refreshes <= sol_refresh.iterations);
    }

    /// If the re-assembly a refresh requires fails, there is no
    /// retreat — the problem has already switched models, so the
    /// system the driver holds describes an objective the problem no
    /// longer evaluates. It must surface, not silently damp a stale
    /// system.
    #[test]
    fn test_refresh_model_reassembly_failure_surfaces() {
        let mut p = AnchoredSurrogate::new(1.0);
        p.next_anchor = Some(3.0);
        p.fail_after_refresh = true;
        let err = solve(&mut p, [0.0], &config(), None).unwrap_err();
        assert!(
            matches!(err, LMError::ModelRefreshFailed { iteration: 1, .. }),
            "{err:?}"
        );
    }

    /// Residual of the model `RefreshOntoStationaryModel` refreshes
    /// ONTO: constant and tiny, so its undamped Gauss-Newton step is
    /// 1e-9 and its quadratic form 1e-18 — stationary to any sane
    /// `qtol` — while no trial can lower its cost.
    const STATIONARY_RESIDUAL: f64 = 1e-9;

    /// A problem that refreshes onto a model which is ALREADY
    /// stationary at the point the refresh happens, but whose trial
    /// costs never improve — so the trial loop can only escalate μ
    /// until the damping cap. This is the shape of a surrogate handing
    /// back to the full model at a converged iterate: the handover
    /// lands on a point the NEW model also calls converged, and every
    /// step it can propose is below the noise.
    struct RefreshOntoStationaryModel {
        refreshed: bool,
    }

    impl RefreshOntoStationaryModel {
        /// Before the refresh, r = x − 1 (so x₀ = 0 is nowhere near
        /// stationary and the iteration-1 battery cannot fire). After
        /// it, the constant stationary residual.
        fn residual(&self, x: f64) -> f64 {
            if self.refreshed {
                STATIONARY_RESIDUAL
            } else {
                x - 1.0
            }
        }
    }

    impl CostProblem<1> for RefreshOntoStationaryModel {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
            let r = self.residual(x[0]);
            Ok(r * r)
        }
        fn refresh_model(&mut self, _x: &[f64; 1]) -> bool {
            if self.refreshed {
                false
            } else {
                self.refreshed = true;
                true
            }
        }
    }

    impl ResidualProblem<1> for RefreshOntoStationaryModel {
        fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
            let r = self.residual(x[0]);
            Ok(NLLSEvaluation {
                residuals: vec![r],
                jacobian: vec![[1.0]],
                cost: r * r,
            })
        }
    }

    /// The convergence battery must re-run against the model a refresh
    /// installed, at the point it was installed at, BEFORE the trial
    /// loop. Without it the driver walks straight into the trials,
    /// exhausts damping and reports FAILURE at a point where the
    /// undamped Gauss-Newton step — the quantity `qtol` is decided on
    /// — is twelve orders of magnitude below the tolerance.
    #[test]
    fn test_refresh_model_convergence_is_tested_on_the_new_model() {
        let mut cfg = config();
        cfg.qtol = 1e-12;
        let mut p = RefreshOntoStationaryModel { refreshed: false };
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();

        assert!(
            sol.converged,
            "declared failure at a point the refreshed model calls converged: {:?}",
            sol.reason
        );
        assert_eq!(sol.reason, TerminationReason::StepTolerance);
        // The refresh happened at x₀ and the driver never left it —
        // the new model was stationary there.
        assert_eq!(sol.x[0], 0.0);
        assert_eq!(sol.iterations, 1);
        // Not one trial was spent: the battery ran first.
        assert_eq!(sol.n_cost_evals, 0);
        // And the reported metric agrees with the reason given.
        let gn = sol.final_gn_qnorm.expect("well-posed after the refresh");
        assert!(gn <= cfg.qtol, "{gn:.3e} vs qtol {:.3e}", cfg.qtol);
    }

    /// The other direction: the post-refresh battery must not fire on
    /// a model that is NOT stationary at the handover point. The
    /// refresh here moves the anchor to 3 while the driver sits at 0,
    /// so the solve must go on to minimize the new objective.
    #[test]
    fn test_refresh_model_does_not_terminate_on_a_moving_model() {
        let mut cfg = config();
        cfg.qtol = 1e-12;
        let mut p = AnchoredSurrogate::new(1.0);
        p.next_anchor = Some(3.0);
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert!((sol.x[0] - 3.0).abs() < 1e-6, "x = {}", sol.x[0]);
        assert!(
            !p.steps("cost").is_empty(),
            "the post-refresh battery swallowed the trial loop"
        );
    }

    /// A problem that accepts steps for one iteration and only THEN
    /// refreshes, onto a model that is stationary where it lands.
    struct RefreshAfterAStep {
        refreshed: bool,
        iterations_seen: usize,
    }

    impl RefreshAfterAStep {
        fn residual(&self, x: f64) -> f64 {
            if self.refreshed {
                STATIONARY_RESIDUAL
            } else {
                x - 100.0
            }
        }
    }

    impl CostProblem<1> for RefreshAfterAStep {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
            let r = self.residual(x[0]);
            Ok(r * r)
        }
        fn refresh_model(&mut self, _x: &[f64; 1]) -> bool {
            self.iterations_seen += 1;
            if self.refreshed || self.iterations_seen < 2 {
                false
            } else {
                self.refreshed = true;
                true
            }
        }
    }

    impl ResidualProblem<1> for RefreshAfterAStep {
        fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
            let r = self.residual(x[0]);
            Ok(NLLSEvaluation {
                residuals: vec![r],
                jacobian: vec![[1.0]],
                cost: r * r,
            })
        }
    }

    /// A model refresh clears the ftol bookkeeping — it must not also
    /// erase the record that a step was ever accepted. The two are
    /// different claims: `actred`/`pred` stop being comparable across a
    /// refresh, but "the last step taken had this size" stays true.
    /// Reporting `None` here says "no step was accepted", which callers
    /// read as "the start point was already stationary" — scott's
    /// `ODResult::update_norm` reports a literal 0.0 on that reading.
    #[test]
    fn test_accepted_step_qnorm_survives_a_model_refresh() {
        let mut cfg = config();
        cfg.qtol = 1e-12;
        let mut p = RefreshAfterAStep {
            refreshed: false,
            iterations_seen: 0,
        };
        let sol = solve(&mut p, [0.0], &cfg, None).unwrap();

        assert!(sol.converged, "{:?}", sol.reason);
        assert!(p.refreshed, "the refresh never happened");
        assert!(
            sol.x[0] != 0.0,
            "the solve must have accepted a step before refreshing"
        );
        let q = sol.accepted_step_qnorm.expect(
            "a solve that accepted a step and then refreshed still accepted a step — \
             None here reads as 'the start point was already stationary'",
        );
        assert!(q > 0.0, "qnorm = {q}");
    }

    // ── Inexact assemblies (served linearizations) ──

    /// A problem that answers the FULL evaluation from a frozen
    /// linearization whenever the driver has just trialled that exact
    /// point, and from its objective otherwise.
    ///
    /// The objective is \\(r(x) = x^2 - 4\\) — genuinely curved, so a
    /// linearization anchored at one point is stationary somewhere the
    /// objective is not, which is the whole hazard. From \\(x_0 = 5\\)
    /// the anchor's model has its root at \\(x = 2.9\\) while the
    /// objective's is at \\(x = 2\\), and the two are far enough apart
    /// that no tolerance can confuse them.
    ///
    /// The serve/don't-serve rule mirrors an orbit fit's: a served
    /// assembly is only ever the REUSE of a trial the solver already
    /// paid for at that point, so the re-assembly the driver forces at a
    /// committed point — where no trial is staged — is exact by
    /// construction.
    struct ServedLinearization {
        /// Where the frozen rows were measured.
        anchor: f64,
        /// The point whose trial cost was last answered from the anchor.
        staged: Option<f64>,
        /// Whether the most recent assembly was served.
        served: bool,
        /// `(x, served)` for every assembly, in order.
        log: Vec<(f64, bool)>,
        /// Report every assembly as exact — the control arm.
        lie: bool,
        /// Coefficient of the truncation allowance a served TRIAL cost
        /// is inflated by, \(K\|x - a\|^2\). An orbit fit carries one
        /// so that a step the linearization cannot vouch for is rejected
        /// rather than taken; large enough, it is what dead-ends a walk.
        allowance: f64,
    }

    impl ServedLinearization {
        fn new() -> Self {
            Self {
                anchor: f64::NAN,
                staged: None,
                served: false,
                log: Vec::new(),
                lie: false,
                allowance: 0.0,
            }
        }
        fn residual(x: f64) -> f64 {
            x * x - 4.0
        }
        fn jacobian(x: f64) -> f64 {
            2.0 * x
        }
        /// The frozen model: \\(r(a) + J(a)\,(x - a)\\).
        fn linearized(&self, x: f64) -> f64 {
            Self::residual(self.anchor) + Self::jacobian(self.anchor) * (x - self.anchor)
        }
        fn served_assemblies(&self) -> usize {
            self.log.iter().filter(|(_, s)| *s).count()
        }
        fn exact_assemblies(&self) -> usize {
            self.log.iter().filter(|(_, s)| !*s).count()
        }
    }

    impl CostProblem<1> for ServedLinearization {
        type Error = TestError;
        fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
            if self.anchor.is_nan() {
                // Nothing frozen yet: the objective answers.
                let r = Self::residual(x[0]);
                return Ok(r * r);
            }
            self.staged = Some(x[0]);
            let r = self.linearized(x[0]);
            let d = x[0] - self.anchor;
            Ok(r * r + self.allowance * d * d)
        }
        fn on_step_accepted(&mut self, _x: &[f64; 1]) {
            self.staged = None;
        }
        fn on_step_rejected(&mut self, _x: &[f64; 1]) {
            self.staged = None;
        }
        fn assembly_is_exact(&self) -> bool {
            self.lie || !self.served
        }
    }

    impl ResidualProblem<1> for ServedLinearization {
        fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
            if self.staged == Some(x[0]) {
                // Reuse what the trial already answered: value from the
                // frozen model, slope the frozen slope.
                self.served = true;
                self.log.push((x[0], true));
                let r = self.linearized(x[0]);
                return Ok(NLLSEvaluation {
                    residuals: vec![r],
                    jacobian: vec![[Self::jacobian(self.anchor)]],
                    cost: r * r,
                });
            }
            self.served = false;
            self.anchor = x[0];
            self.log.push((x[0], false));
            let r = Self::residual(x[0]);
            Ok(NLLSEvaluation {
                residuals: vec![r],
                jacobian: vec![[Self::jacobian(x[0])]],
                cost: r * r,
            })
        }
    }

    /// THE guarantee: a convergence verdict may only latch on an exact
    /// assembly, so the point a solve certifies is the OBJECTIVE's
    /// stationary point and not a linearization's.
    ///
    /// The control arm is the same problem lying about the same
    /// assemblies, and it lands 0.9 away — at the frozen model's own
    /// root. That gap is what the hook buys; without it the driver
    /// cannot tell the two apart, because on the served system every
    /// convergence criterion is genuinely satisfied.
    #[test]
    fn test_convergence_never_latches_on_an_inexact_assembly() {
        let mut honest = ServedLinearization::new();
        let sol = solve(&mut honest, [5.0], &config(), None).unwrap();

        assert!(sol.converged, "{:?}", sol.reason);
        assert!(
            (sol.x[0] - 2.0).abs() < 1e-6,
            "converged at {} — the objective's root is 2",
            sol.x[0]
        );
        assert!(
            honest.served_assemblies() > 0,
            "nothing was served, so the guard was never exercised"
        );
        assert!(
            !honest.log.last().expect("assemblies were logged").1,
            "the last assembly before the verdict was served"
        );
        assert!(
            sol.n_model_refreshes > 0,
            "no re-assembly was forced, so no walk was ever ended"
        );

        let mut lying = ServedLinearization::new();
        lying.lie = true;
        let fiction = solve(&mut lying, [5.0], &config(), None).unwrap();
        assert!(fiction.converged, "{:?}", fiction.reason);
        assert!(
            (fiction.x[0] - 2.9).abs() < 1e-6,
            "a problem that vouches for its linearization should land at ITS root (2.9), \
             got {}",
            fiction.x[0]
        );
    }

    /// The default is `true`, so a problem that never mentions the hook
    /// sees the driver it always saw: one assembly per accepted point,
    /// no forced re-assemblies, and a verdict off the system it built.
    #[test]
    fn test_assembly_is_exact_defaults_to_true_and_forces_nothing() {
        let mut p = Tracked::new(overshoot_rational);
        let sol = solve(&mut p, [6.5], &config(), None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        assert_eq!(sol.n_model_refreshes, 0);
    }

    /// An exhaustion reached on a served system is a fact about the
    /// SERVED system. With the truncation allowance turned up, the walk
    /// across one anchor can only end that way: a couple of steps in,
    /// the allowance outgrows the reduction the frozen model can still
    /// promise, every further trial is rejected, and the driver runs out
    /// of damping — all while the served residual is nowhere near zero,
    /// so no convergence criterion can fire on it either.
    ///
    /// The driver must re-anchor and carry on. Reporting
    /// `DampingExhausted` here would be reporting the stand-in's dead end
    /// as the objective's, at a point 0.9 away from the objective's root.
    #[test]
    fn test_an_exhausted_stand_in_re_anchors_instead_of_ending_the_solve() {
        let mut p = ServedLinearization::new();
        p.allowance = 50.0;
        // One accepted step per anchor is all the allowance permits, so
        // the walk contracts by ~0.75 an iteration and needs a budget
        // three times the default to certify. That IS the quasi-Newton
        // trade this hook exists to make legible: cheap steps, more of
        // them.
        let mut cfg = config();
        cfg.max_iterations = 300;
        let sol = solve(&mut p, [5.0], &cfg, None).unwrap();

        assert!(sol.converged, "{:?}", sol.reason);
        assert!(
            (sol.x[0] - 2.0).abs() < 1e-6,
            "the objective's root is 2, the solve stopped at {}",
            sol.x[0]
        );
        assert!(
            p.served_assemblies() > 0 && p.exact_assemblies() > 1,
            "served {} exact {} — the walk-then-re-anchor pattern never happened",
            p.served_assemblies(),
            p.exact_assemblies()
        );
        assert!(
            sol.n_rejected_trials > 0,
            "no trial was ever rejected, so no walk dead-ended"
        );
        assert!(sol.n_model_refreshes >= 2, "{}", sol.n_model_refreshes);
    }

    /// The iteration budget is the one stop the driver does not choose,
    /// so it is the one that can land mid-walk — and every [`LMSolution`]
    /// field a caller reads off the final system would then describe the
    /// stand-in. That matters most where it is least visible: a consumer
    /// scoring its stopping point in units of the fit's own covariance
    /// (scott's converged-in-practice stall acceptance does exactly this)
    /// would be dividing the linearization's q-norm by the
    /// linearization's covariance, which is self-consistent and therefore
    /// always small. Any walk endpoint would rubber-stamp.
    ///
    /// So the final re-assembly is not optional either. For
    /// \\(r(x) = x^2 - 4\\) the exact system's undamped Gauss-Newton
    /// q-norm has the closed form \\((x^2-4)^2\\), which the frozen model
    /// at any other anchor does not reproduce — so this pins the number
    /// and not merely the code path.
    #[test]
    fn test_the_returned_diagnostics_are_exact_even_out_of_budget() {
        let mut p = ServedLinearization::new();
        p.allowance = 50.0;
        let mut cfg = config();
        // Small enough that the budget runs out while a walk is in
        // progress, which the 300-iteration arm above is chosen to avoid.
        cfg.max_iterations = 12;
        let sol = solve(&mut p, [5.0], &cfg, None).unwrap();

        assert!(!sol.converged, "{:?}", sol.reason);
        assert!(
            !p.log.last().expect("assemblies were logged").1,
            "the solve returned on a served assembly"
        );
        let x = sol.x[0];
        let want = (x * x - 4.0) * (x * x - 4.0);
        let got = sol.final_gn_qnorm.expect("a non-singular 1x1 system");
        assert!(
            (got - want).abs() <= 1e-12 * want.max(1e-12),
            "final_gn_qnorm {got} is not the objective's {want} at x = {x}"
        );
    }

    /// The forced re-assembly is not optional. A problem that serves a
    /// second linearization there leaves the driver with no system it
    /// may judge and no honest way to stop, and that must surface.
    #[test]
    fn test_a_permanently_inexact_assembly_surfaces() {
        struct AlwaysServes;
        impl CostProblem<1> for AlwaysServes {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok((x[0] - 2.0) * (x[0] - 2.0))
            }
            fn assembly_is_exact(&self) -> bool {
                false
            }
        }
        impl ResidualProblem<1> for AlwaysServes {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 2.0],
                    jacobian: vec![[1.0]],
                    cost: (x[0] - 2.0) * (x[0] - 2.0),
                })
            }
        }
        let err = solve(&mut AlwaysServes, [10.0], &config(), None).unwrap_err();
        assert!(
            matches!(err, LMError::InexactAssemblyPersisted { .. }),
            "{err:?}"
        );
    }

    // ── Reported step norms ──

    /// `accepted_step_qnorm` and `final_gn_qnorm` measure different
    /// things and only the latter is comparable to `qtol`. On a
    /// μ-starved solve the accepted (damped) step is vanishingly
    /// small at a garbage iterate while the undamped Gauss-Newton step
    /// — the quantity the test is actually decided on — is enormous.
    /// Quoting the former against a tolerance is what makes a failure
    /// message read as "converged".
    #[test]
    fn test_final_gn_qnorm_is_the_convergence_metric_not_the_accepted_step() {
        fn stiff_cost(x: f64) -> f64 {
            if x.abs() < 1e-7 {
                1e6 - 1e-3 * (1e-7 - x.abs()) / 1e-7
            } else {
                1e6 + 1.0
            }
        }
        struct StiffValley;
        impl CostProblem<1> for StiffValley {
            type Error = TestError;
            fn evaluate_cost(&mut self, x: &[f64; 1]) -> Result<f64, TestError> {
                Ok(stiff_cost(x[0]))
            }
        }
        impl ResidualProblem<1> for StiffValley {
            fn evaluate(&mut self, x: &[f64; 1]) -> Result<NLLSEvaluation<1>, TestError> {
                Ok(NLLSEvaluation {
                    residuals: vec![x[0] - 1000.0],
                    jacobian: vec![[1.0]],
                    cost: stiff_cost(x[0]),
                })
            }
        }
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.xtol = 0.0;
        cfg.ftol = 0.0;
        cfg.qtol = 1.0;
        cfg.max_iterations = 30;
        let sol = solve(&mut StiffValley, [0.0], &cfg, None).unwrap();
        assert!(!sol.converged, "{:?}", sol.reason);

        let gn = sol.final_gn_qnorm.expect("undamped system is well posed");
        assert!(
            gn > cfg.qtol,
            "the metric convergence is decided on must be ABOVE the tolerance on a \
             non-converged solve: {gn:.3e} vs qtol {:.3e}",
            cfg.qtol
        );
        if let Some(accepted) = sol.accepted_step_qnorm {
            assert!(
                accepted < cfg.qtol,
                "this regression needs a μ-crushed accepted step to be meaningful: {accepted:.3e}"
            );
            assert!(
                gn / accepted > 1e6,
                "the two norms should be orders of magnitude apart: {gn:.3e} vs {accepted:.3e}"
            );
        }
    }

    /// On a converged solve the reported metric is below the tolerance
    /// it was tested against — the other direction of the same
    /// contract.
    #[test]
    fn test_final_gn_qnorm_below_qtol_on_a_step_tolerance_exit() {
        let mut cfg = config();
        cfg.gtol = 0.0;
        cfg.ftol = 0.0;
        cfg.qtol = 1e-12;
        let mut p = Tracked::new(overshoot_rational);
        let sol = solve(&mut p, [6.5], &cfg, None).unwrap();
        assert!(sol.converged, "{:?}", sol.reason);
        let gn = sol.final_gn_qnorm.expect("well-posed at the solution");
        assert!(gn <= cfg.qtol, "{gn:.3e} vs qtol {:.3e}", cfg.qtol);
    }

    const GOLDEN_MU_FINAL_BITS: u64 = 4565320297239836560;
    const GOLDEN_ACCEPTED_BITS: [u64; 10] = [
        4619004367821864960, // x0 = 6.5 (the initial commit)
        4618792513637290406,
        4618291282969069438,
        4617770008120803146,
        4617507755895746747,
        4617364428888646038,
        4617321199543857098,
        4617315762137427473,
        4617315521566082796,
        4617315517979513644, // = GOLDEN_X_BITS (the returned x)
    ];
    const GOLDEN_REJECTED_BITS: [u64; 5] = [
        4594659349414467680,
        4594887588787221536,
        4596247515049879712,
        4602933352259637776,
        4614750054244392936,
    ];

    // x = 5.0000000159096025, cost ≈ 6.3e-15 — the overshoot_rational
    // problem from x0 = 6.5 under the default config, traced on macOS
    // arm64. Pure fixed-order f64 arithmetic end to end, so these bits
    // are the contract on every IEEE 754 platform.
    const GOLDEN_X_BITS: u64 = 4617315517979513644;
    const GOLDEN_COST_BITS: u64 = 4394527585940473101;
    const GOLDEN_ITERATIONS: usize = 10;
    const GOLDEN_COST_EVALS: usize = 14;
    const GOLDEN_REJECTED: usize = 5;
}