spectrafit_solver/

dispatch.rs

//! `fit()` — the public solver entry point and dispatcher.
//!
//! Parses [`FitOptionsSpec::solver`] into a [`Solver`] and routes it: `irls`,
//! `global` and `varpro` dispatch early to their own module/crate, while the
//! Levenberg–Marquardt family (`lm`/`trf`/`geodesic`/`lm-legacy`, and `auto`'s
//! non-VarPro case) runs the shared pre-solve, solve and post-fit path here.
//!
//! Steps (LM-family path):
//!   1. Compile the graph → get `free_keys`
//!   2. Extract initial parameter vector + bounds
//!   3. Evaluate the model *before* optimisation (`init_fit`)
//!   4. Run `LevenbergMarquardt::minimize`
//!   5. Compute post-fit statistics: χ², reduced-χ², R², DOF, AIC, BIC
//!   6. Compute covariance matrix: (J_wᵀ J_w)⁻¹ when σ is user-supplied, else (JᵀJ)⁻¹ · (χ²/DOF)
//!   7. Assemble and return `FitResultSpec`

use std::cell::RefCell;
use std::collections::HashMap;
use std::ops::ControlFlow;
use std::str::FromStr;

use levenberg_marquardt::LevenbergMarquardt;
use nalgebra::DVector;
use spectrafit_graph::{evaluate_compiled, CompiledGraph};
use spectrafit_types::{
    CoreError, FitGraphSpec, FitOptionsSpec, FitResultSpec, MeasurementSpec, TerminationReason,
};

use crate::error::SolverError;
use crate::global::{solve_global, DeConfig};
use crate::irls::{solve_irls, WeightFn, DEFAULT_MAX_OUTER_ITER, DEFAULT_TOL_WEIGHTS};
use crate::lm_problem::LmProblem;
use crate::postfit::{self, LmSolveOutcome};

/// Flatten one dataset's `dims × points` coordinates into the **point-major**
/// (stride = n_dims) layout the executor's `coord_layout` expects: point `i`
/// occupies `x[i*n_dims .. (i+1)*n_dims]`. For the 1-D case this is exactly
/// `ds.x[0]`, so existing behaviour is unchanged; for n-D it interleaves the
/// per-dimension coordinate rows so a Gaussian2D node receives `[xᵢ, yᵢ]`.
pub(crate) fn point_major_x(ds: &MeasurementSpec) -> Vec<f64> {
    let n_dims = ds.x.len();
    if n_dims == 0 {
        return Vec::new();
    }
    let n_points = ds.x[0].len();
    let mut out = Vec::with_capacity(n_points * n_dims);
    for i in 0..n_points {
        for dim in &ds.x {
            out.push(dim[i]);
        }
    }
    out
}

// ---------------------------------------------------------------------------
// Solver dispatch
// ---------------------------------------------------------------------------

/// Every solver [`fit`] can dispatch to, parsed once from
/// [`FitOptionsSpec::solver`]. This enum + [`Solver::parse`] + the single match
/// in [`fit`] is the *one place* that answers "which solvers exist and where
/// does each one run":
///
/// * [`Lm`](Solver::Lm) / [`Trf`](Solver::Trf) / [`Geodesic`](Solver::Geodesic) /
///   [`LmLegacy`](Solver::LmLegacy) — the Levenberg–Marquardt family. All run the
///   shared pre-solve + solve below: the first three on the faer
///   `spectrafit-levenberg-marquardt` core (TRF = Coleman–Li bound scaling,
///   geodesic = Transtrum acceleration), `LmLegacy` on the nalgebra oracle.
/// * [`Irls`](Solver::Irls) / [`Global`](Solver::Global) / [`Varpro`](Solver::Varpro)
///   — distinct algorithms that dispatch early to their own module/crate.
/// * [`Auto`](Solver::Auto) — VarPro when the graph is separable, unconstrained
///   and untied (see [`graph_prefers_varpro`]), otherwise the faer LM family.
#[derive(Debug, Clone, Copy, PartialEq)]
enum Solver {
    /// Default faer Levenberg–Marquardt.
    Lm,
    /// nalgebra `levenberg-marquardt` parity oracle.
    LmLegacy,
    /// Trust-Region-Reflective (Coleman–Li bound scaling on the faer LM core).
    Trf,
    /// Geodesic acceleration on the faer LM core.
    Geodesic,
    /// Powell's dogleg trust-region method.
    Dogleg,
    /// Matrix-free Newton-CG (Steihaug–Toint) trust-region method.
    NewtonCg,
    /// Iteratively-reweighted least squares (robust loss).
    Irls(WeightFn),
    /// Differential Evolution (global search).
    Global,
    /// Variable Projection (separable nonlinear least squares).
    Varpro,
    /// Choose from the graph shape at run time.
    Auto,
}

impl Solver {
    /// Parse the `FitOptionsSpec.solver` string.
    ///
    /// Unrecognised strings are rejected rather than silently defaulting to
    /// LM — previously a typo (e.g. `"lmm"`, wrong-case `"Trf"`) silently ran
    /// the LM solver instead of erroring, which could mask a misconfigured
    /// fit with a plausible-looking but unintended result.
    fn parse(s: &str) -> Result<Self, SolverError> {
        match s {
            "lm" => Ok(Solver::Lm),
            "lm-legacy" => Ok(Solver::LmLegacy),
            "trf" => Ok(Solver::Trf),
            "geodesic" | "lm-geodesic" => Ok(Solver::Geodesic),
            "dogleg" => Ok(Solver::Dogleg),
            "newton-cg" | "newton_cg" | "newtoncg" | "steihaug" => Ok(Solver::NewtonCg),
            "global" => Ok(Solver::Global),
            "varpro" => Ok(Solver::Varpro),
            "auto" => Ok(Solver::Auto),
            _ if s.starts_with("irls") => {
                let name = s.split_once(':').map(|(_, name)| name).unwrap_or("huber");
                Ok(Solver::Irls(WeightFn::from_str(name)?))
            }
            other => Err(SolverError::UnrecognisedSolver(other.to_string())),
        }
    }
}

/// Whether a graph carries any tied parameters — a constraint expression from a
/// graph `expr_edge` OR a per-parameter `Parameter.expr`. VarPro reconstructs
/// parameters from the linear/nonlinear split and never applies the `TiedPlan`,
/// so a tied graph must never be routed to / accepted by VarPro on EITHER surface
/// (CX-VPE-01: checking `expr_edges` alone silently dropped `Parameter.expr` ties).
fn graph_has_tied_params(graph: &FitGraphSpec) -> bool {
    !graph.expr_edges.is_empty()
        || graph
            .nodes
            .iter()
            .any(|n| n.parameters.values().any(|p| p.expr.is_some()))
}

/// Whether `solver="auto"` should choose VarPro: the graph is separable, every
/// nonlinear parameter is unconstrained (VarPro ignores bounds), and there are no
/// tied parameters (via `expr_edges` or `Parameter.expr`; VarPro cannot honour them).
fn graph_prefers_varpro(graph: &FitGraphSpec) -> bool {
    !graph_has_tied_params(graph)
        // VarPro stacks all datasets and is not dataset_index-aware, so it must
        // not be auto-selected for simultaneous multi-dataset ("global analysis")
        // graphs — those need the LM-family scoped executor path.
        && graph.nodes.iter().all(|n| n.dataset_index.is_none())
        && spectrafit_varpro::is_separable(graph)
        && graph.nodes.iter().all(|n| {
            n.parameters
                .iter()
                // amplitude is the linear coeff VarPro projects out; skip it BY NAME
                // (HashMap order is non-deterministic, so positional `.skip(1)` could
                // drop an arbitrary nonlinear param and miss its bounds — silent bound
                // violation when that param is bounded but amplitude is not).
                .filter(|(name, _)| name.as_str() != "amplitude")
                .filter(|(_, p)| p.vary)
                .all(|(_, p)| p.min.is_infinite() && p.max.is_infinite())
        })
}

/// Build the `ParameterSpec` map and run VarPro. Shared by the explicit
/// `solver="varpro"` arm and the `solver="auto"`→VarPro arm.
fn solve_varpro_path(
    graph: &FitGraphSpec,
    datasets: &[MeasurementSpec],
    options: &FitOptionsSpec,
) -> Result<FitResultSpec, CoreError> {
    let mut param_specs: HashMap<String, spectrafit_types::ParameterSpec> = HashMap::new();
    for node in &graph.nodes {
        for (pname, pspec) in &node.parameters {
            param_specs.insert(format!("{}.{}", node.id, pname), pspec.clone());
        }
    }
    spectrafit_varpro::solve_varpro(graph, datasets, &param_specs, options)
}

// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------

/// Run a fit of `graph` against `datasets` using the solver named by
/// `options.solver` (see `Solver::parse` for the accepted strings: `lm`,
/// `lm-legacy`, `trf`, `geodesic`, `dogleg`, `newton-cg`, `irls[:weight]`,
/// `global`, `varpro`, `auto`).
///
/// The LM family (`lm`, `lm-legacy`, `trf`, `geodesic`) plus the
/// trust-region methods (`dogleg`, `newton-cg`) and `auto`'s non-VarPro
/// fallback share one pre-solve/solve/post-fit path through this function.
/// `irls`, `global`, and `varpro` are routed to their own modules/crates in
/// `dispatch_solver` and return before that shared path runs.
///
/// Returns a fully populated [`FitResultSpec`] including parameter values,
/// standard errors, covariance matrix, and goodness-of-fit statistics.
///
/// # Errors
/// Returns [`CoreError`] when:
/// - `options.solver` does not parse (unrecognised solver name, or an
///   unrecognised `irls:<weight>` weight-function name) — `Solver::parse`.
/// - The routed non-LM-family solver fails: `irls` (`crate::irls::solve_irls`),
///   `global` (differential evolution), or `varpro` (not separable, tied
///   parameters present, or per-dataset scoping unsupported).
/// - Graph compilation or the initial model evaluation fails (malformed
///   graph, a free-key/parameter lookup miss, or a model evaluation error) —
///   `prepare_lm_pre_solve`.
/// - The LM-family solve itself fails or produces a numerically unusable
///   result (e.g. ill-conditioned Jacobian) that maps to `CoreError::Solver`.
/// - `options.eta` is set alongside `solver="dogleg"` / `"newton-cg"` but is
///   not finite and in `[0.0, 0.25)` — `SolverError::InvalidEta`.
/// - Post-fit assembly fails while re-evaluating the model or its Jacobian
///   at the solution — `postfit::assemble_result`.
pub fn fit(
    graph: &FitGraphSpec,
    datasets: Vec<MeasurementSpec>,
    options: &FitOptionsSpec,
) -> Result<FitResultSpec, CoreError> {
    // Run faer's dense kernels serially for the whole fit (solve loop AND the
    // post-fit covariance/condition-number). Our matrices are skinny (p ≲ 150),
    // so faer's default Rayon dispatch costs more in thread-pool wake-up
    // (~30–50 µs per matmul) than the arithmetic itself — measured ~4–6×
    // slowdown overall, ~4.7× on a 50k×3 fit, including the n×p post-fit
    // Gram/SVD. The data-parallelism that pays off lives one level down in the
    // graph executor's own Rayon path (over points), not in these
    // per-iteration p×p factorizations.
    //
    // THIS IS THE ONLY SITE that sets faer's global parallelism. The LM and
    // trust-region drivers each used to set it again on entry to their own
    // `minimize`/`minimize_tr`; setting one process-global three times from
    // three crates made the effective policy depend on call order and left the
    // drivers quietly mutating global state their doc comments never mentioned.
    // Ownership now sits here, at the one place that already owns the whole
    // fit — faer is solver-internal in this workspace, so a single global
    // policy set at the fit boundary is safe.
    faer::set_global_parallelism(faer::Par::Seq);

    // ── 0. Solver dispatch ───────────────────────────────────────────────────
    // Route the non-LM-family solvers to their own module/crate; the LM family
    // (Lm/LmLegacy/Trf/Geodesic, and Auto's non-VarPro case) falls through to the
    // shared pre-solve + solve below.

    let solver = Solver::parse(&options.solver)?;
    let datasets = match dispatch_solver(solver, graph, datasets, options) {
        ControlFlow::Break(result) => return result,
        ControlFlow::Continue(datasets) => datasets,
    };

    // ── 1/2/3/6/7. Compile the graph, collect free-param specs, evaluate the
    // initial model, and build the index-based node-param buffers ───────────
    let LmPreSolve {
        cg,
        free_keys,
        all_params,
        bounds,
        scales,
        sigma,
        x_all,
        y_all,
        init_fit,
        init_params,
        node_param_bufs,
        free_to_node_param,
        tied_to_node_param,
    } = prepare_lm_pre_solve(graph, &datasets)?;

    // ── 8. Construct problem ─────────────────────────────────────────────────
    let problem = LmProblem {
        compiled: &cg,
        free_keys: free_keys.clone(),
        bounds: bounds.clone(),
        all_params: all_params.clone(),
        node_param_bufs,
        free_to_node_param,
        tied_to_node_param,
        x_concat: x_all.clone(),
        y_concat: y_all.clone(),
        params: init_params,
        scales: scales.clone(),
        sigma,
        residual_buf: RefCell::new(vec![0.0; y_all.len()]),
        // One Jacobian row per data point (residual), not per coordinate: size by
        // the point count (y_all.len()), which equals x_all.len() only in 1-D.
        jacobian_buf: RefCell::new(vec![0.0; y_all.len() * free_keys.len()]),
        residual_count: RefCell::new(0),
        jacobian_count: RefCell::new(0),
        residual_time_ns: RefCell::new(0),
        jacobian_time_ns: RefCell::new(0),
    };

    // ── 4/8. Initialise tied parameters and run the LM-family solve ──────────
    let outcome = run_lm_solve(problem, solver, options, free_keys.len(), y_all.len())?;

    // ── 5/6/7. Post-fit statistics, covariance, and result assembly ──────────
    postfit::assemble_result(
        outcome,
        postfit::PostfitInputs {
            cg: &cg,
            graph,
            datasets: &datasets,
            x_all: &x_all,
            y_all: &y_all,
        },
        init_fit,
    )
}

// ---------------------------------------------------------------------------
// Private helpers
// ---------------------------------------------------------------------------

/// Route `solver` to its own module/crate when it is not part of the LM
/// family, otherwise hand `datasets` back unchanged. Step 0 of [`fit`]'s
/// sequence, extracted verbatim from the original inline `match`.
///
/// Returns [`ControlFlow::Break`] carrying the already-final
/// `Result<FitResultSpec, CoreError>` for `Irls` / `Global` / `Varpro` / an
/// `Auto` that prefers VarPro — `fit` returns this immediately. Returns
/// [`ControlFlow::Continue`] carrying `datasets` handed back unchanged for the
/// LM family (`Lm` / `LmLegacy` / `Trf` / `Geodesic` / `Dogleg` / `NewtonCg`,
/// and an `Auto` that does not prefer VarPro) — `fit` falls through to the
/// shared pre-solve + solve path for these.
fn dispatch_solver(
    solver: Solver,
    graph: &FitGraphSpec,
    datasets: Vec<MeasurementSpec>,
    options: &FitOptionsSpec,
) -> ControlFlow<Result<FitResultSpec, CoreError>, Vec<MeasurementSpec>> {
    match solver {
        Solver::Irls(weight) => ControlFlow::Break(solve_irls(
            graph,
            datasets,
            options,
            weight,
            DEFAULT_MAX_OUTER_ITER,
            // The IRLS weight-convergence tolerance is its own constant, not
            // `options.tolerance` (the LM cost/step/gradient tolerance) — see
            // `DEFAULT_TOL_WEIGHTS`.
            DEFAULT_TOL_WEIGHTS,
        )),

        Solver::Global => {
            // Honour the caller's iteration budget: map `max_iterations` onto the
            // number of DE generations. Without this the 500-generation budget the
            // benchmark passes for pathological cases was silently dropped and DE
            // ran with the 100-generation default.
            let config = DeConfig {
                max_gen: if options.max_iterations > 0 {
                    options.max_iterations as usize
                } else {
                    DeConfig::default().max_gen
                },
                ..DeConfig::default()
            };
            ControlFlow::Break(solve_global(graph, datasets, options, config))
        }
        Solver::Varpro => {
            // Explicit VarPro: require separability, and reject tied params from
            // EITHER surface (expr_edges or Parameter.expr) — VarPro reconstructs
            // parameters from the linear/nonlinear split and never applies the
            // tied-plan, so a tie would be silently dropped (CX-VPE-01).
            if !spectrafit_varpro::is_separable(graph) {
                return ControlFlow::Break(Err(SolverError::VarproNotSeparable.into()));
            }
            if graph_has_tied_params(graph) {
                return ControlFlow::Break(Err(SolverError::VarproExprEdgesUnsupported.into()));
            }
            // VarPro stacks all datasets into one synthetic spec and is not
            // dataset_index-aware, so a per-dataset (local) node would silently
            // contribute to every dataset. Reject rather than mis-scope; the
            // LM-family solvers honour dataset_index (simultaneous global analysis).
            if graph.nodes.iter().any(|n| n.dataset_index.is_some()) {
                return ControlFlow::Break(
                    Err(SolverError::VarproDatasetScopingUnsupported.into()),
                );
            }
            // Falls through to the shared VarPro tail below, which applies the
            // post-fit success guards. Do not call `solve_varpro_path` directly
            // here — see the note on that tail.
            ControlFlow::Break(run_varpro_guarded(graph, &datasets, options))
        }
        Solver::Auto if graph_prefers_varpro(graph) => {
            // MUST use the same guarded tail as the explicit `Solver::Varpro`
            // arm above. This arm previously called `solve_varpro_path` raw,
            // skipping `finalize_varpro_result` — so on a graph where automatic
            // routing chooses VarPro, `solver="auto"` reported success on
            // exactly the degenerate fits `solver="varpro"` reported as failed,
            // on identical data. `auto` is the documented default, which made
            // it a silent wrong answer on the recommended path.
            ControlFlow::Break(run_varpro_guarded(graph, &datasets, options))
        }
        // LM family + the Δ-radius methods (and Auto without VarPro) flow through
        // to the shared pre-solve + solve below.
        Solver::Lm
        | Solver::LmLegacy
        | Solver::Trf
        | Solver::Geodesic
        | Solver::Dogleg
        | Solver::NewtonCg
        | Solver::Auto => ControlFlow::Continue(datasets),
    }
}

/// Bundle of everything the shared LM-family pre-solve step produces (steps
/// 1/2/3/6/7 of [`fit`]'s sequence: compile the graph, collect free-parameter
/// specs, evaluate the initial model, and build the index-based node-param
/// buffers). [`fit`] destructures this into the same flat local bindings the
/// original inline code used, so field names mirror the original variable
/// names exactly.
struct LmPreSolve {
    /// Compiled DAG (step 1), including `dataset_offsets` for simultaneous
    /// multi-dataset ("global analysis") fits.
    cg: CompiledGraph,
    /// `"node_id.param_name"` for every free parameter, canonical order.
    free_keys: Vec<String>,
    /// All parameters (free + fixed), keyed by `"node_id.param_name"`.
    all_params: HashMap<String, f64>,
    /// `(min, max)` bounds per free parameter.
    bounds: Vec<(f64, f64)>,
    /// Per-free-parameter `Parameter.scale` factors (identity when unset).
    scales: Vec<f64>,
    /// Per-point σ (default 1.0 where `MeasurementSpec.sigma` is `None`).
    sigma: Vec<f64>,
    /// Point-major concatenated x across all datasets.
    x_all: Vec<f64>,
    /// Concatenated y across all datasets.
    y_all: Vec<f64>,
    /// The model evaluated at the initial parameters, before optimisation.
    init_fit: Vec<f64>,
    /// Initial free-parameter vector in the optimiser's scaled working units.
    init_params: DVector<f64>,
    /// Per-node parameter buffers (fixed pre-filled, free updated in-place).
    node_param_bufs: Vec<Vec<f64>>,
    /// free-param index → `(node_idx, param_pos)`.
    free_to_node_param: Vec<(usize, usize)>,
    /// tied-target index → `(node_idx, param_pos)`.
    tied_to_node_param: Vec<(usize, usize)>,
}

/// Compile `graph`, collect free-parameter specs, evaluate the initial model,
/// and build the index-based node-param buffers — steps 1/2/3/6/7 of [`fit`]'s
/// sequence, extracted verbatim from the original inline code (including the
/// `Parameter.scale` wiring doc comment and its `debug_assert_eq!`).
fn prepare_lm_pre_solve(
    graph: &FitGraphSpec,
    datasets: &[MeasurementSpec],
) -> Result<LmPreSolve, CoreError> {
    // ── 1. Compile graph ────────────────────────────────────────────────────
    let mut cg = CompiledGraph::compile(graph)?;
    // Per-dataset point boundaries (cumulative, len = n_datasets + 1) so the
    // executor can scope a node with `dataset_index = Some(i)` to dataset i's
    // contiguous point-range — the primitive for simultaneous multi-dataset
    // ("global analysis") fits. Left as the single-range no-op for one dataset.
    cg.dataset_offsets = {
        let mut offs = Vec::with_capacity(datasets.len() + 1);
        let mut acc = 0usize;
        offs.push(0);
        for ds in datasets {
            acc += ds.y.len();
            offs.push(acc);
        }
        offs
    };
    let free_keys = cg.free_keys.clone();

    // ── 2/3. All-param flat map + initial values/bounds/scales for free params ─
    let (all_params, init_vals, bounds, scales) = collect_free_param_specs(graph, &free_keys)?;

    let init_vals_clone = init_vals.clone();
    // The solver optimises the scaled working variable θ'_i = θ_i / s_i, so its
    // starting point is the physical initial value divided by the scale factor
    // (a no-op when s_i = 1.0). `node_param_bufs` is still seeded from the
    // physical `init_vals_clone` below.
    let init_params = DVector::from_vec(
        init_vals
            .iter()
            .zip(scales.iter())
            .map(|(&v, &s)| v / s)
            .collect::<Vec<f64>>(),
    );

    // ── 4. Build per-point sigma vector (default = 1.0) ─────────────────────
    let sigma: Vec<f64> = datasets
        .iter()
        .flat_map(|ds| {
            let n = ds.y.len();
            match &ds.sigma {
                Some(s) => s.clone(),
                None => vec![1.0_f64; n],
            }
        })
        .collect();

    // ── 5. Concatenated x and y for all datasets (cached in the problem) ────
    // Point-major (stride = n_dims) so the executor reads each point's full
    // coordinate; 1-D reduces to the old per-dataset `ds.x[0]` concatenation.
    let x_all: Vec<f64> = datasets.iter().flat_map(point_major_x).collect();

    let y_all: Vec<f64> = datasets
        .iter()
        .flat_map(|ds| ds.y.iter().copied())
        .collect();

    // ── 6. Evaluate initial model (before optimisation) ─────────────────────
    let init_fit = evaluate_compiled(&cg, &all_params, &x_all)?;

    // ── 7. Build index-based node param buffers ──────────────────────────────
    let (node_param_bufs, free_to_node_param, tied_to_node_param) =
        build_node_param_buffers(&cg, &all_params, &bounds, &init_vals_clone)?;

    // ── 7b. Per-parameter scale normalization (Parameter.scale → LM step) ─────
    //
    // `Parameter.scale` lets a caller request that a parameter be optimised in
    // rescaled units (e.g. an amplitude ~1e6 alongside a width ~1e-3). The solver
    // operates on the rescaled variable `θ'_i = θ_i / s_i`, which equalises the
    // columns of the Jacobian and thereby reshapes — and typically improves — the
    // conditioning of `JᵀJ` (see the condition-number computation in §11). The
    // wiring is now end-to-end:
    //   1. `init_params` (above) is the physical start divided by `scales[i]`,
    //   2. `LmProblem` carries `scales`; `set_params`/`params` translate between
    //      scaled (optimiser) and physical (`node_param_bufs`) coordinates and the
    //      analytic/FD Jacobian columns are multiplied by `scales[i]` (chain rule),
    //   3. `node_param_bufs` — and hence `to_flat` / the reported result — is
    //      always physical, so the converged values come back in real units.
    // The post-fit condition number (§11 of postfit) applies the same column
    // scaling so it reflects the optimiser's effective conditioning.
    //
    // With all-default (None ⇒ 1.0) scales every step above is an exact
    // arithmetic no-op, so a fit without any `scale` set is byte-for-byte
    // unchanged — the parity oracle (`tests/parity.rs`) pins this.
    debug_assert_eq!(
        scales.len(),
        free_keys.len(),
        "one scale factor per free parameter"
    );

    Ok(LmPreSolve {
        cg,
        free_keys,
        all_params,
        bounds,
        scales,
        sigma,
        x_all,
        y_all,
        init_fit,
        init_params,
        node_param_bufs,
        free_to_node_param,
        tied_to_node_param,
    })
}

/// Initialise tied parameters from the starting free values, then dispatch
/// `solver` to the matching LM-family engine — extracted verbatim from the
/// original inline code: `lm-legacy` on the retained nalgebra
/// `levenberg-marquardt` crate (the parity oracle); `dogleg` / `newton-cg` on
/// the Δ-radius trust-region drivers; everything else (`lm` / `trf` /
/// `geodesic` / `auto`) on the faer-native LM core.
fn run_lm_solve<'a>(
    mut problem: LmProblem<'a>,
    solver: Solver,
    options: &FitOptionsSpec,
    free_keys_len: usize,
    y_all_len: usize,
) -> Result<LmSolveOutcome<'a>, CoreError> {
    // Initialise tied parameters from the starting free values. `node_param_bufs`
    // was seeded with the tied-target *placeholders* (their spec `value`, e.g.
    // 0.0); without this the first residual/Jacobian would evaluate the model
    // with stale tied values. No-op for graphs without ties (`expr_edges` or `Parameter.expr`).
    if problem.has_tied() {
        let p0: Vec<f64> = problem.params.iter().copied().collect();
        problem.set_free_and_tied(&p0);
    }

    // `patience` controls max function evaluations = patience * (n_params + 1)
    let patience = (options.max_iterations as usize).max(1);

    // Solver engine selection. The faer-native trust-region core is the default
    // for every LM-family solver string (`"lm"`, `"auto"`, `"trf"`, `"geodesic"`);
    // the proven `levenberg-marquardt` (nalgebra) crate is retained as the
    // `"lm-legacy"` oracle for parity testing until it is removed.
    let tol = if options.tolerance > 0.0 {
        options.tolerance
    } else {
        1e-8
    };
    let max_nfev = patience.saturating_mul(free_keys_len + 1);
    let _solve_t0 = std::time::Instant::now();
    let (
        result_problem,
        n_iter_val,
        success_val,
        message_val,
        cost_history,
        gradient_norm_history,
        params_history,
    ) = match solver {
        Solver::LmLegacy => {
            // Retained nalgebra `levenberg-marquardt` crate — the parity oracle.
            // It does not expose a per-iteration trajectory, so the history stays
            // empty (the benchmark layer reconstructs a labelled proxy).
            let lm = LevenbergMarquardt::new().with_patience(patience);
            let lm = if options.tolerance > 0.0 {
                lm.with_tol(options.tolerance)
            } else {
                lm
            };
            let (rp, report) = lm.minimize(problem);
            (
                rp,
                report.number_of_evaluations as u64,
                report.termination.was_successful(),
                map_termination(&report.termination).as_str().to_string(),
                Vec::new(),
                Vec::new(),
                Vec::new(),
            )
        }
        Solver::Dogleg | Solver::NewtonCg => {
            // Explicit Δ-radius trust-region methods on the shared framework loop.
            // `Report`/`Termination` are the same types the faer LM core returns.
            // Power-user knobs for the TR core: each `Option<f64>` on
            // `FitOptionsSpec` keeps the library default when `None`, so the wire
            // surface remains backward-compatible. The TR driver clamps `delta0`
            // and `max_delta` internally; we just pass through what the caller set.
            let mut cfg = spectrafit_dogleg::TrustRegionConfig {
                ftol: tol,
                xtol: tol,
                gtol: tol,
                max_nfev,
                ..Default::default()
            };
            if let Some(d0) = options.delta0 {
                cfg.delta0 = d0;
            }
            if let Some(md) = options.max_delta {
                cfg.max_delta = md;
            }
            if let Some(e) = options.eta {
                // `TrustRegionConfig::eta` is only `debug_assert!`-checked by
                // the driver, which compiles out in release builds — this is
                // the check that actually enforces the bound for callers that
                // reach the Rust core directly (the JSON `fit` entrypoint),
                // not just the Pydantic `FitOptions` validator.
                if !(e.is_finite() && (0.0..0.25).contains(&e)) {
                    return Err(SolverError::InvalidEta(e).into());
                }
                cfg.eta = e;
            }
            let report = if solver == Solver::Dogleg {
                spectrafit_dogleg::minimize(&mut problem, &cfg)
            } else {
                spectrafit_newton_cg::minimize(&mut problem, &cfg)
            };
            let msg = faer_termination_str(report.termination).to_string();
            (
                problem,
                report.n_iter as u64,
                report.termination.was_successful(),
                msg,
                report.cost_history,
                report.gradient_norm_history,
                report.params_history,
            )
        }
        // LM family: lm / trf / geodesic / auto on the faer LM core. Regime-adaptive
        // (normal equations for tall-skinny, SVD/secular for many-parameter), with
        // Moré column scaling applied inside the driver.
        _ => {
            let cfg = spectrafit_levenberg_marquardt::StrategyConfig {
                kind: spectrafit_levenberg_marquardt::select_regime(y_all_len, free_keys_len),
                ftol: tol,
                xtol: tol,
                gtol: tol,
                max_nfev,
                // `solver="geodesic"` (Transtrum acceleration) — faster on sloppy /
                // degenerate multi-peak surfaces; otherwise plain faer LM.
                geodesic: solver == Solver::Geodesic,
                // `solver="trf"` (Trust-Region-Reflective) — Coleman–Li bound scaling
                // so steps shrink near active bounds. `solver="auto"` that did not
                // route to VarPro also lands here as TRF: across a representative
                // case per scenario family, TRF is the fastest LM-family strategy at
                // top accuracy in most problem classes. The TRF-over-VarPro speed
                // result may invert for very large separable problems, which is why
                // VarPro routing is retained for the fully-unconstrained separable
                // case.
                bound_scaling: solver == Solver::Trf || solver == Solver::Auto,
                ..Default::default()
            };
            let report = spectrafit_levenberg_marquardt::minimize(&mut problem, &cfg);
            let msg = faer_termination_str(report.termination).to_string();
            (
                problem,
                report.n_iter as u64,
                report.termination.was_successful(),
                msg,
                report.cost_history,
                report.gradient_norm_history,
                report.params_history,
            )
        }
    };
    let solve_ns = _solve_t0.elapsed().as_nanos();

    Ok(LmSolveOutcome {
        result_problem,
        n_iter_val,
        success_val,
        message_val,
        solve_ns,
        cost_history,
        gradient_norm_history,
        params_history,
    })
}

/// `(all_params, init_vals, bounds, scales)` — see [`collect_free_param_specs`].
type FreeParamSpecs = (HashMap<String, f64>, Vec<f64>, Vec<(f64, f64)>, Vec<f64>);

/// `(node_param_bufs, free_to_node_param, tied_to_node_param)` — see
/// [`build_node_param_buffers`].
type NodeParamBuffers = (Vec<Vec<f64>>, Vec<(usize, usize)>, Vec<(usize, usize)>);

/// Build the flat all-parameters map (free + fixed) plus per-free-parameter
/// initial values, bounds, and scale factors.
///
/// Returns `(all_params, init_vals, bounds, scales)`, all indexed in `free_keys`
/// order for the latter three. `scales[i]` defaults to 1.0 (identity — a
/// non-finite or non-positive `Parameter.scale` is meaningless for
/// normalization and falls back to the identity too).
fn collect_free_param_specs(
    graph: &FitGraphSpec,
    free_keys: &[String],
) -> Result<FreeParamSpecs, CoreError> {
    let mut all_params: HashMap<String, f64> = HashMap::new();
    for node in &graph.nodes {
        for (pname, pspec) in &node.parameters {
            let key = format!("{}.{}", node.id, pname);
            all_params.insert(key, pspec.value);
        }
    }

    let mut init_vals: Vec<f64> = Vec::with_capacity(free_keys.len());
    let mut bounds: Vec<(f64, f64)> = Vec::with_capacity(free_keys.len());
    let mut scales: Vec<f64> = Vec::with_capacity(free_keys.len());

    for key in free_keys {
        // key = "node_id.param_name"
        let (node_id, param_name) = key
            .split_once('.')
            .ok_or_else(|| SolverError::Dispatch(format!("malformed free key: '{}'", key)))?;

        let node = graph
            .nodes
            .iter()
            .find(|n| n.id == node_id)
            .ok_or_else(|| SolverError::Dispatch(format!("node '{}' not found", node_id)))?;

        let pspec = node.parameters.get(param_name).ok_or_else(|| {
            SolverError::Dispatch(format!(
                "param '{}' not found in node '{}'",
                param_name, node_id
            ))
        })?;

        init_vals.push(pspec.value);
        bounds.push((pspec.min, pspec.max));
        // A non-finite or non-positive scale is meaningless for normalization;
        // fall back to the identity (1.0) so it is a no-op.
        let s = pspec.scale.unwrap_or(1.0);
        scales.push(if s.is_finite() && s > 0.0 { s } else { 1.0 });
    }

    Ok((all_params, init_vals, bounds, scales))
}

/// Build index-based node param buffers, the free→node/param index mapping,
/// and the tied-target→node/param index mapping.
///
/// `node_param_bufs[i]` holds `cg.nodes[i]`'s current param values in model
/// `param_names()` order — fixed params are pre-filled and never changed;
/// free params are updated in-place by `LmProblem::set_params()`.
/// `free_to_node_param` inverts `cg.node_free_cols` (built during `compile()`)
/// so free-key index → `(node_idx, param_pos_in_node)`, and is used here to
/// apply the initial bounds-clamp to `node_param_bufs`. `tied_to_node_param`
/// maps each tied target (`tied_plan.order`) to its `(node_idx, param_pos)`
/// slot so the solver can write each recomputed tied value back — empty when
/// the graph declares no ties (`expr_edges` or `Parameter.expr`).
fn build_node_param_buffers(
    cg: &CompiledGraph,
    all_params: &HashMap<String, f64>,
    bounds: &[(f64, f64)],
    init_vals_clone: &[f64],
) -> Result<NodeParamBuffers, CoreError> {
    let mut node_param_bufs: Vec<Vec<f64>> = (0..cg.nodes.len())
        .map(|i| {
            cg.node_params(i, all_params)
                .unwrap_or_else(|_| vec![0.0; cg.nodes[i].param_names.len()])
        })
        .collect();

    let free_to_node_param: Vec<(usize, usize)> = {
        // Invert node_free_cols: for each free_key index (col), find the
        // (node_idx, local_param_idx) pair.
        let mut mapping = vec![(0usize, 0usize); cg.free_keys.len()];
        for (node_idx, pairs) in cg.node_free_cols.iter().enumerate() {
            for &(local_idx, col) in pairs {
                mapping[col] = (node_idx, local_idx);
            }
        }
        mapping
    };

    // Apply initial bounds clamping to node_param_bufs for free params.
    for (i, &(node_idx, param_pos)) in free_to_node_param.iter().enumerate() {
        let (lo, hi) = bounds[i];
        node_param_bufs[node_idx][param_pos] = init_vals_clone[i].clamp(lo, hi);
    }

    let tied_to_node_param: Vec<(usize, usize)> = cg
        .tied_plan
        .order
        .iter()
        .map(|tp| {
            let (nid, pname) = tp
                .target
                .split_once('.')
                .ok_or_else(|| SolverError::MalformedTiedTarget(tp.target.clone()))?;
            let ni = cg
                .nodes
                .iter()
                .position(|n| n.id == nid)
                .ok_or_else(|| SolverError::TiedTargetNodeMissing(nid.to_string()))?;
            let pos = cg.nodes[ni]
                .param_names
                .iter()
                .position(|p| p == pname)
                .ok_or_else(|| SolverError::TiedTargetParamMissing(pname.to_string()))?;
            Ok::<(usize, usize), CoreError>((ni, pos))
        })
        .collect::<Result<_, _>>()?;

    Ok((node_param_bufs, free_to_node_param, tied_to_node_param))
}

/// Run the VarPro path and apply the shared post-fit success guards.
///
/// Both VarPro entry points — the explicit `Solver::Varpro` arm and the
/// `Solver::Auto` arm when `graph_prefers_varpro` — must route through here.
/// They are separate match arms because they have different *preconditions*
/// (explicit VarPro rejects non-separable/tied/dataset-scoped graphs with a
/// specific error; automatic routing only reaches VarPro when the graph already
/// qualifies), but they must share the same *postcondition*. Keeping the guard
/// in one function is what makes "the two solver strings agree on `success`"
/// structural rather than a thing to remember.
fn run_varpro_guarded(
    graph: &FitGraphSpec,
    datasets: &[MeasurementSpec],
    options: &FitOptionsSpec,
) -> Result<FitResultSpec, CoreError> {
    let result = solve_varpro_path(graph, datasets, options)?;
    Ok(finalize_varpro_result(graph, datasets, result))
}

/// Apply the shared post-fit success guards (off-domain runaway + degenerate
/// peak-collapse) to a VarPro result.
///
/// VarPro builds its own [`FitResultSpec`] outside `postfit::assemble_result`
/// (it has no `LmProblem`/covariance path to hang off), so this re-derives the
/// inputs `apply_postfit_guards` needs — free-key list, flat parameter map,
/// concatenated x/y — from the already-assembled `result` and applies the same
/// guard the LM-family path gets for free. Without this a degenerate VarPro
/// fit would report false success.
fn finalize_varpro_result(
    graph: &FitGraphSpec,
    datasets: &[MeasurementSpec],
    mut result: FitResultSpec,
) -> FitResultSpec {
    let x_all: Vec<f64> = datasets.iter().flat_map(point_major_x).collect();
    let y_all: Vec<f64> = datasets
        .iter()
        .flat_map(|ds| ds.y.iter().copied())
        .collect();
    let final_flat: HashMap<String, f64> = result
        .parameters
        .iter()
        .map(|(k, p)| (k.clone(), p.value))
        .collect();
    let vp_free_keys: Vec<String> = result
        .parameters
        .iter()
        .filter(|(_, p)| p.vary)
        .map(|(k, _)| k.clone())
        .collect();
    let (success, message) = postfit::apply_postfit_guards(
        graph,
        &vp_free_keys,
        &final_flat,
        &x_all,
        &y_all,
        result.r_squared,
        vp_free_keys.len(),
        result.success,
        result.message.clone(),
    );
    result.success = success;
    result.message = message;
    result
}

/// Stable snake_case message for a faer-native [`spectrafit_levenberg_marquardt::Termination`].
fn faer_termination_str(t: spectrafit_levenberg_marquardt::Termination) -> &'static str {
    use spectrafit_levenberg_marquardt::Termination as T;
    match t {
        T::Gtol => "converged_gtol",
        T::Ftol => "converged_ftol",
        T::Xtol => "converged_xtol",
        T::ResidualsZero => "residuals_zero",
        T::MaxEval => "max_iterations",
        T::NoImprovement => "no_improvement_possible",
        T::NumericalError => "numerical_error",
    }
}

/// Map the `levenberg_marquardt` crate's `TerminationReason` onto our own
/// stable enum so callers get a deterministic snake_case string, not a
/// debug-format blob that changes with crate versions.
fn map_termination(lm: &levenberg_marquardt::TerminationReason) -> TerminationReason {
    use levenberg_marquardt::TerminationReason as L;
    match lm {
        L::ResidualsZero => TerminationReason::ResidualsZero,
        L::Orthogonal => TerminationReason::Orthogonal,
        L::Converged { .. } => TerminationReason::Converged,
        L::LostPatience => TerminationReason::MaxIterations,
        L::NoImprovementPossible(_) => TerminationReason::NoImprovementPossible,
        L::NoParameters => TerminationReason::NoParameters,
        L::NoResiduals => TerminationReason::NoResiduals,
        L::WrongDimensions(_) => TerminationReason::WrongDimensions,
        L::Numerical(_) => TerminationReason::NumericalError,
        L::User(_) => TerminationReason::UserCancelled,
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------
#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use spectrafit_types::{
        FitGraphSpec, FitOptionsSpec, MeasurementSpec, ModelNodeSpec, ModelTypeStr, ParameterSpec,
    };
    use std::collections::HashMap;

    // ── helpers ──────────────────────────────────────────────────────────────

    fn make_param(value: f64, vary: bool) -> ParameterSpec {
        ParameterSpec {
            value,
            min: f64::NEG_INFINITY,
            max: f64::INFINITY,
            vary,
            expr: None,
            scale: None,
        }
    }

    fn default_options() -> FitOptionsSpec {
        FitOptionsSpec {
            schema_version: None,
            solver: "lm".to_string(),
            max_iterations: 200,
            tolerance: 1e-8,
            delta0: None,
            max_delta: None,
            eta: None,
        }
    }

    /// Analytical Gaussian: A · exp(−(x−c)² / (2σ²))
    fn gaussian(x: f64, amplitude: f64, center: f64, sigma: f64) -> f64 {
        amplitude * (-(x - center).powi(2) / (2.0 * sigma * sigma)).exp()
    }

    // ── Test 1: Gaussian parameter recovery ──────────────────────────────────

    #[test]
    fn test_gaussian_recovery() {
        // True parameters
        let (true_a, true_c, true_s) = (5.0_f64, 2.0_f64, 0.5_f64);

        // x grid: 50 points in [−1, 5]
        let n = 50usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -1.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        let y: Vec<f64> = x
            .iter()
            .map(|&xi| gaussian(xi, true_a, true_c, true_s))
            .collect();

        // Graph with perturbed initial params
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(4.0, true));
        params.insert("center".into(), make_param(1.8, true));
        params.insert("sigma".into(), make_param(0.6, true));

        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };

        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };

        let result = fit(&graph, vec![dataset], &default_options()).expect("fit should not error");

        assert!(
            result.success,
            "LM should converge; message: {}",
            result.message
        );
        assert!(
            result.n_iter > 0,
            "should have done at least one evaluation"
        );
        assert!(
            result.chi2 < 1e-10,
            "chi2 = {} should be near zero",
            result.chi2
        );

        let a = result.parameters["g1.amplitude"].value;
        let c = result.parameters["g1.center"].value;
        let s = result.parameters["g1.sigma"].value;

        assert_relative_eq!(a, true_a, max_relative = 0.01);
        assert_relative_eq!(c, true_c, max_relative = 0.01);
        assert_relative_eq!(s, true_s, max_relative = 0.01);
    }

    // ── Test: unrecognised `solver`/weight-fn strings are rejected, not
    //    silently defaulted ────────────────────────────────────────────────
    //
    // Before this fix, `Solver::parse` silently mapped ANY unrecognised
    // string to `Solver::Lm` (and `WeightFn::from_str` silently mapped any
    // unrecognised name to Huber) — a typo like `solver="lmm"` ran a fit
    // with the LM solver instead of erroring, which could hide a
    // misconfigured request behind a plausible-looking successful result.

    #[test]
    fn fit_rejects_unrecognised_solver_string_instead_of_defaulting_to_lm() {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(4.0, true));
        params.insert("center".into(), make_param(1.8, true));
        params.insert("sigma".into(), make_param(0.6, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![vec![0.0, 1.0, 2.0]],
            y: vec![1.0, 2.0, 1.0],
            sigma: None,
            label: None,
        };
        let opts = FitOptionsSpec {
            solver: "lmm".to_string(), // typo for "lm"
            ..default_options()
        };

        let err = fit(&graph, vec![dataset], &opts)
            .expect_err("a typo'd solver name must error, not silently run LM");
        let msg = err.to_string();
        assert!(
            msg.contains("lmm"),
            "error should name the offending string, got: {msg}"
        );
        assert!(
            msg.contains("expected one of"),
            "error should list valid solver names, got: {msg}"
        );
    }

    // ── Test: `eta >= 0.25` is rejected on the JSON/Rust path, not just by
    //    the Pydantic `Field(ge=0.0, lt=0.25)` validator ─────────────────────
    //
    // Before this fix, `FitOptionsSpec.eta` was only bounds-checked in
    // `python/spectrafit_core/options.py`; a caller reaching the Rust core
    // directly (the JSON `fit` entrypoint, or any non-Python consumer) could
    // set `eta >= 0.25`, which the trust-region driver only caught with a
    // `debug_assert!` that compiles out of release builds — so a release
    // build silently ran with the invalid ratio.

    #[test]
    fn fit_rejects_eta_at_or_above_a_quarter_on_dogleg() {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(4.0, true));
        params.insert("center".into(), make_param(1.8, true));
        params.insert("sigma".into(), make_param(0.6, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![vec![0.0, 1.0, 2.0]],
            y: vec![1.0, 2.0, 1.0],
            sigma: None,
            label: None,
        };
        let opts = FitOptionsSpec {
            solver: "dogleg".to_string(),
            eta: Some(0.5),
            ..default_options()
        };

        let err = fit(&graph, vec![dataset], &opts)
            .expect_err("eta >= 0.25 must error, not silently run with an invalid ratio");
        let msg = err.to_string();
        assert!(
            msg.contains("eta"),
            "error should name the offending option, got: {msg}"
        );
    }

    #[test]
    fn fit_rejects_unrecognised_irls_weight_fn_instead_of_defaulting_to_huber() {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(4.0, true));
        params.insert("center".into(), make_param(1.8, true));
        params.insert("sigma".into(), make_param(0.6, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![vec![0.0, 1.0, 2.0]],
            y: vec![1.0, 2.0, 1.0],
            sigma: None,
            label: None,
        };
        let opts = FitOptionsSpec {
            solver: "irls:buisquare".to_string(), // typo for "bisquare"
            ..default_options()
        };

        let err = fit(&graph, vec![dataset], &opts)
            .expect_err("a typo'd irls weight-fn name must error, not silently run Huber");
        let msg = err.to_string();
        assert!(
            msg.contains("buisquare"),
            "error should name the offending string, got: {msg}"
        );
    }

    // ── Test 2: Constant model trivial fit ────────────────────────────────────

    #[test]
    fn test_constant_recovery() {
        let n = 20usize;
        let x: Vec<f64> = (0..n).map(|i| i as f64 / (n - 1) as f64).collect();
        // Tiny deterministic "noise" to avoid a perfectly flat problem
        let y: Vec<f64> = x
            .iter()
            .enumerate()
            .map(|(i, _)| 3.0 + 1e-6 * (i as f64 * 0.1).sin())
            .collect();

        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("c".into(), make_param(0.0, true)); // start far from 3

        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "const1".into(),
                model_type: ModelTypeStr::Constant,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };

        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };

        let result = fit(&graph, vec![dataset], &default_options()).expect("fit should not error");

        assert!(
            result.success,
            "LM should converge; message: {}",
            result.message
        );

        let c_val = result.parameters["const1.c"].value;
        assert!(
            (c_val - 3.0).abs() < 1e-4,
            "recovered constant = {}, expected ≈ 3.0",
            c_val
        );
    }

    // ── SP-2: N-D (gaussian_nd) end-to-end fits ───────────────────────────────
    //
    // The compiler infers D from the node's `center_<i>` parameters and builds a
    // `GaussianND` of that dimensionality; the N-D-general executor strides the
    // coordinate buffer at stride D and the LM solver recovers every parameter.
    // Both 3-D and 5-D run, so "arbitrary N" is exercised, not extrapolated.
    fn run_gaussian_nd_recovery(d: usize) {
        let amp = 3.0_f64;
        let centers: Vec<f64> = (0..d).map(|i| -1.0 + 0.5 * i as f64).collect();
        let sigmas: Vec<f64> = (0..d).map(|i| 0.8 + 0.1 * i as f64).collect();

        // Keep the grid tiny so d=5 stays cheap (5^5 = 3125 points).
        let axis_n = if d <= 3 { 7 } else { 5 };
        let axis: Vec<f64> = (0..axis_n)
            .map(|i| -3.0 + 6.0 * i as f64 / (axis_n - 1) as f64)
            .collect();
        // Cartesian product of the axis over d dimensions → point coordinates.
        let mut coords: Vec<Vec<f64>> = vec![vec![]];
        for _dim in 0..d {
            let mut next = Vec::with_capacity(coords.len() * axis.len());
            for prefix in &coords {
                for &a in &axis {
                    let mut p = prefix.clone();
                    p.push(a);
                    next.push(p);
                }
            }
            coords = next;
        }
        let g = |pt: &[f64]| -> f64 {
            let mut z = 0.0;
            for i in 0..d {
                let dx = pt[i] - centers[i];
                z -= dx * dx / (2.0 * sigmas[i] * sigmas[i]);
            }
            amp * z.exp()
        };
        let y: Vec<f64> = coords.iter().map(|pt| g(pt)).collect();
        let n_points = coords.len();
        // Dimension-major x: x[dim] is the coordinate array for axis `dim`.
        let x: Vec<Vec<f64>> = (0..d)
            .map(|dim| (0..n_points).map(|p| coords[p][dim]).collect())
            .collect();

        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(2.0, true));
        for (i, &center) in centers.iter().enumerate() {
            params.insert(format!("center_{i}"), make_param(center + 0.3, true));
            params.insert(format!("sigma_{i}"), make_param(1.0, true));
        }
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "gnd".into(),
                model_type: ModelTypeStr::GaussianNd,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let dataset = MeasurementSpec {
            schema_version: None,
            x,
            y,
            sigma: None,
            label: None,
        };
        let result =
            fit(&graph, vec![dataset], &default_options()).expect("N-D fit should not error");
        assert!(
            result.success,
            "LM should converge on {d}-D; msg: {}",
            result.message
        );
        let p = &result.parameters;
        assert_relative_eq!(p["gnd.amplitude"].value, amp, max_relative = 0.02);
        for i in 0..d {
            assert_relative_eq!(
                p[&format!("gnd.center_{i}")].value,
                centers[i],
                epsilon = 0.05
            );
            assert_relative_eq!(
                p[&format!("gnd.sigma_{i}")].value,
                sigmas[i],
                epsilon = 0.05
            );
        }
    }

    #[test]
    fn gaussian_nd_fit_recovers_3d() {
        run_gaussian_nd_recovery(3);
    }

    #[test]
    fn gaussian_nd_fit_recovers_5d_arbitrary_n() {
        run_gaussian_nd_recovery(5);
    }

    // ── TDD: tied-parameter (expr_edge) end-to-end fit ────────────────────────
    //
    // The compiler builds a dependency-ordered `tied_plan` (parse +
    // topo-order + cycle-detection). The LM solver loop calls `TiedPlan::apply`
    // on every iteration, so tied parameters are recomputed from `expr_edges` and
    // `Parameter.expr` at each step. These tests run unconditionally and pass.

    /// Build a two-Gaussian graph where `g2.amplitude = k * g1.amplitude`.
    ///
    /// The tie is expressed through the top-level `expr_edge` list only.
    /// `ParameterSpec.expr` is intentionally `None` on the tied amplitude so
    /// this helper exercises the legacy `expr_edge` path without triggering the
    /// `DuplicateExprTarget` error that T1 introduced when both routes target
    /// the same parameter.
    fn tied_two_gaussian_graph(k: f64) -> FitGraphSpec {
        let mut g1: HashMap<String, ParameterSpec> = HashMap::new();
        g1.insert("amplitude".into(), make_param(4.0, true));
        g1.insert("center".into(), make_param(-1.0, true));
        g1.insert("sigma".into(), make_param(0.5, true));

        let mut g2: HashMap<String, ParameterSpec> = HashMap::new();
        // `vary = false`, `expr = None` — the tie lives in expr_edges below.
        g2.insert("amplitude".into(), make_param(0.0, false));
        g2.insert("center".into(), make_param(1.0, true));
        g2.insert("sigma".into(), make_param(0.5, true));

        FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![
                ModelNodeSpec {
                    id: "g1".into(),
                    model_type: ModelTypeStr::Gaussian,
                    dataset_index: None,
                    parameters: g1,
                },
                ModelNodeSpec {
                    id: "g2".into(),
                    model_type: ModelTypeStr::Gaussian,
                    dataset_index: None,
                    parameters: g2,
                },
            ],
            expr_edges: vec![spectrafit_types::ExprEdge {
                target_node: "g2".into(),
                target_param: "amplitude".into(),
                expression: format!("{} * g1.amplitude", k),
            }],
        }
    }

    /// A tied fit recovers `g2.amplitude == k * g1.amplitude` (wired in M6).
    #[test]
    fn test_tied_amplitude_fit_recovers_ratio() {
        let k = 0.5_f64;
        let (true_a, true_s) = (5.0_f64, 0.5_f64);
        let n = 80usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -3.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        let y: Vec<f64> = x
            .iter()
            .map(|&xi| gaussian(xi, true_a, -1.0, true_s) + gaussian(xi, k * true_a, 1.0, true_s))
            .collect();

        let graph = tied_two_gaussian_graph(k);
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        let result = fit(&graph, vec![dataset], &default_options()).unwrap();

        assert!(
            result.success,
            "tied fit should converge: {}",
            result.message
        );
        let a1 = result.parameters["g1.amplitude"].value;
        let a2 = result.parameters["g2.amplitude"].value;
        // The tied plan enforces the ratio exactly each iteration…
        assert_relative_eq!(a2, k * a1, epsilon = 1e-9);
        // …and the free amplitude is recovered.
        assert_relative_eq!(a1, true_a, max_relative = 1e-3);
    }

    /// End-to-end proof that a `Parameter.expr`-only tie (no `expr_edge`) is
    /// honoured by the solver — the tied parameter tracks its derived value at
    /// the converged solution, not its initial placeholder.
    ///
    /// The graph expresses `g2.amplitude = 0.5 * g1.amplitude` purely through
    /// `ParameterSpec.expr`; the `expr_edges` list is intentionally **empty**.
    /// T1 folds `Parameter.expr` into the compiled `TiedPlan`, and the solver
    /// calls `set_free_and_tied` each iteration, so `g2.amplitude` is recomputed
    /// from the current `g1.amplitude` on every step.
    #[test]
    fn test_param_expr_fit_recovers_derived_value() {
        let k = 0.5_f64;
        let (true_a, true_s) = (5.0_f64, 0.5_f64);
        let n = 80usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -3.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        // Synthesise data from the known-true parameters consistent with the tie.
        let y: Vec<f64> = x
            .iter()
            .map(|&xi| gaussian(xi, true_a, -1.0, true_s) + gaussian(xi, k * true_a, 1.0, true_s))
            .collect();

        // Build a graph identical to `tied_two_gaussian_graph(k)` **except** the
        // `expr_edges` list is empty — the tie lives solely in `Parameter.expr`.
        let mut g1: HashMap<String, ParameterSpec> = HashMap::new();
        g1.insert("amplitude".into(), make_param(4.0, true));
        g1.insert("center".into(), make_param(-1.0, true));
        g1.insert("sigma".into(), make_param(0.5, true));

        let mut g2: HashMap<String, ParameterSpec> = HashMap::new();
        let mut tied_amp = make_param(0.0, false);
        tied_amp.expr = Some(format!("{} * g1.amplitude", k));
        g2.insert("amplitude".into(), tied_amp);
        g2.insert("center".into(), make_param(1.0, true));
        g2.insert("sigma".into(), make_param(0.5, true));

        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![
                ModelNodeSpec {
                    id: "g1".into(),
                    model_type: ModelTypeStr::Gaussian,
                    dataset_index: None,
                    parameters: g1,
                },
                ModelNodeSpec {
                    id: "g2".into(),
                    model_type: ModelTypeStr::Gaussian,
                    dataset_index: None,
                    parameters: g2,
                },
            ],
            // Intentionally empty — the tie is expressed only via Parameter.expr.
            expr_edges: vec![],
        };

        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        let result = fit(&graph, vec![dataset], &default_options()).unwrap();

        assert!(
            result.success,
            "param-expr tied fit should converge: {}",
            result.message
        );
        let a1 = result.parameters["g1.amplitude"].value;
        let a2 = result.parameters["g2.amplitude"].value;
        // `g2.amplitude` must equal `k * g1.amplitude` at the solution — not
        // the initial placeholder (0.0).
        assert!(
            a2.abs() > 1e-6,
            "tied param must not be frozen at placeholder 0.0"
        );
        assert_relative_eq!(a2, k * a1, epsilon = 1e-9);
        // The free amplitude must recover the ground truth.
        assert_relative_eq!(a1, true_a, max_relative = 1e-3);
    }

    /// The tied parameter is excluded from the free set, so the solver reports
    /// DOF = n_points − n_free with n_free reduced by the tied count.
    #[test]
    fn test_tied_fit_reduces_free_param_count() {
        let graph = tied_two_gaussian_graph(0.5);
        let cg = CompiledGraph::compile(&graph).unwrap();
        // 6 params total, 1 tied (g2.amplitude) → 5 free.
        assert_eq!(cg.free_keys.len(), 5);
        assert_eq!(cg.tied_plan.len(), 1);
    }

    // ── Helper: build a well-conditioned single-Gaussian fit problem ──────────

    /// Returns a graph + dataset that recovers a clean Gaussian. `scale` is a
    /// base `Parameter.scale` factor (`None` ⇒ all parameters unset, i.e. 1.0).
    ///
    /// When `Some(base)`, a deliberately **non-uniform** scale is applied across
    /// the three parameters — `amplitude·base`, `center` unscaled, `sigma/base`.
    /// This is intentional: a *uniform* column scaling `J → J·(s·I)` multiplies
    /// `JᵀJ` by `s²` and leaves `κ(JᵀJ) = σ_max/σ_min` exactly invariant, so it
    /// could never change the reported conditioning. Spreading the scale across
    /// the columns changes their relative norms, which is what actually reshapes
    /// `κ(JᵀJ)` (see `parameter_scale_changes_effective_conditioning`).
    fn gaussian_fit_inputs(scale: Option<f64>) -> (FitGraphSpec, MeasurementSpec) {
        let (true_a, true_c, true_s) = (5.0_f64, 2.0_f64, 0.5_f64);
        let n = 50usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -1.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        let y: Vec<f64> = x
            .iter()
            .map(|&xi| gaussian(xi, true_a, true_c, true_s))
            .collect();

        let mk = |value: f64, scale: Option<f64>| ParameterSpec {
            value,
            min: f64::NEG_INFINITY,
            max: f64::INFINITY,
            vary: true,
            expr: None,
            scale,
        };
        // Non-uniform scale spread (see doc comment) so κ(JᵀJ) genuinely changes.
        let (sa, sc, ss) = match scale {
            None => (None, None, None),
            Some(base) => (Some(base), Some(1.0), Some(1.0 / base)),
        };
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), mk(4.0, sa));
        params.insert("center".into(), mk(1.8, sc));
        params.insert("sigma".into(), mk(0.6, ss));

        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        (graph, dataset)
    }

    // ── GREEN: condition number is computed end-to-end through the LM path ────

    #[test]
    fn condition_number_is_some_and_finite_after_gaussian_fit() {
        let (graph, dataset) = gaussian_fit_inputs(None);
        let result = fit(&graph, vec![dataset], &default_options()).expect("fit should not error");
        let cond = result
            .condition_number
            .expect("well-conditioned Gaussian fit should report a condition number");
        assert!(cond.is_finite(), "condition number must be finite: {cond}");
        assert!(cond >= 1.0, "condition number must be ≥ 1.0: {cond}");
    }

    // ── U3: Parameter.scale is wired into the LM step / Jacobian ──────────────

    #[test]
    fn parameter_scale_changes_effective_conditioning() {
        // Same problem, once unscaled and once with a deliberately non-uniform
        // `Parameter.scale` spread across the free parameters (see
        // `gaussian_fit_inputs`). With scale wired into the LM step and Jacobian,
        // the reported condition number must differ — the per-column scaling
        // reshapes the columns of J and hence κ(JᵀJ). (A *uniform* scale would be
        // κ-invariant, which is why the helper spreads it across columns.)
        let (g_unscaled, d_unscaled) = gaussian_fit_inputs(None);
        let (g_scaled, d_scaled) = gaussian_fit_inputs(Some(1000.0));

        let r_unscaled =
            fit(&g_unscaled, vec![d_unscaled], &default_options()).expect("unscaled fit");
        let r_scaled = fit(&g_scaled, vec![d_scaled], &default_options()).expect("scaled fit");

        let c0 = r_unscaled
            .condition_number
            .expect("unscaled fit should report κ");
        let c1 = r_scaled
            .condition_number
            .expect("scaled fit should report κ");
        assert!(
            (c0 - c1).abs() > 1e-6,
            "Parameter.scale should change effective conditioning: κ_unscaled={c0}, κ_scaled={c1}"
        );
    }

    #[test]
    fn parameter_scale_of_one_is_a_bitwise_no_op() {
        // The parity contract: `scale = Some(1.0)` on every free parameter must be
        // byte-for-byte identical to `scale = None` (the un-scaled path). Every
        // scaling operation reduces to multiply/divide by 1.0, which is exact in
        // IEEE-754, so the reported κ — and the recovered parameters — must match
        // to the last bit, not merely within a tolerance.
        let (g_none, d_none) = gaussian_fit_inputs(None);
        let g_one = {
            let mut g = g_none.clone();
            for p in g.nodes[0].parameters.values_mut() {
                p.scale = Some(1.0);
            }
            g
        };

        let r_none = fit(&g_none, vec![d_none.clone()], &default_options()).expect("none fit");
        let r_one = fit(&g_one, vec![d_none], &default_options()).expect("scale=1 fit");

        assert_eq!(
            r_none.condition_number, r_one.condition_number,
            "scale=Some(1.0) must reproduce the un-scaled κ bit-for-bit"
        );
        for (key, p_none) in &r_none.parameters {
            let p_one = &r_one.parameters[key];
            assert_eq!(
                p_none.value, p_one.value,
                "scale=1 changed converged value of {key}: {} vs {}",
                p_none.value, p_one.value
            );
        }
    }

    // ── auto-routing: non-VarPro graph → TRF ─────────────────────────────────

    #[test]
    fn auto_routes_to_trf_for_bounded_graph() {
        // A bounded sigma disqualifies the graph from VarPro (graph_prefers_varpro
        // requires unconstrained nonlinear params), so `solver="auto"` must fall
        // through to TRF. Assert auto's result is identical to an explicit trf fit
        // (and that graph_prefers_varpro is indeed false for this graph).
        let (true_a, true_c, true_s) = (5.0_f64, 2.0_f64, 0.5_f64);
        let n = 60usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -1.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        let y: Vec<f64> = x
            .iter()
            .map(|&xi| gaussian(xi, true_a, true_c, true_s))
            .collect();

        let bounded = |value: f64, min: f64, max: f64| ParameterSpec {
            value,
            min,
            max,
            vary: true,
            expr: None,
            scale: None,
        };
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), bounded(4.0, 0.0, f64::INFINITY));
        params.insert("center".into(), make_param(1.8, true));
        params.insert("sigma".into(), bounded(0.6, 1e-6, 10.0)); // finite → not varpro-eligible
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        assert!(
            !graph_prefers_varpro(&graph),
            "bounded sigma must disqualify the graph from VarPro"
        );
        let data = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        let opts = |s: &str| FitOptionsSpec {
            solver: s.to_string(),
            ..default_options()
        };
        let r_auto = fit(&graph, vec![data.clone()], &opts("auto")).expect("auto fit");
        let r_trf = fit(&graph, vec![data], &opts("trf")).expect("trf fit");
        // auto routed to TRF → bit-for-bit the same outcome.
        assert_eq!(
            r_auto.n_iter, r_trf.n_iter,
            "auto should match trf iterations"
        );
        assert_relative_eq!(r_auto.chi2, r_trf.chi2, max_relative = 1e-12);
        for key in ["g1.amplitude", "g1.center", "g1.sigma"] {
            assert_relative_eq!(
                r_auto.parameters[key].value,
                r_trf.parameters[key].value,
                max_relative = 1e-12
            );
        }
    }

    // ── degenerate peak-collapse guard ───────────────────────────────────────

    #[test]
    fn degenerate_peak_collapse_is_flagged_unsuccessful() {
        // Narrow Gaussian (true centre 7) fit from centre=0 on the flat tail:
        // local LM stalls and collapses the amplitude to ~0 (R² < 0). The
        // degenerate-fit guard must downgrade success to false.
        let n = 200usize;
        let x: Vec<f64> = (0..n).map(|i| 10.0 * i as f64 / (n - 1) as f64).collect();
        let y: Vec<f64> = x.iter().map(|&xi| gaussian(xi, 3.0, 7.0, 0.3)).collect();

        let bounded = |value: f64, min: f64, max: f64| ParameterSpec {
            value,
            min,
            max,
            vary: true,
            expr: None,
            scale: None,
        };
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), bounded(1.0, 0.0, f64::INFINITY));
        params.insert("center".into(), make_param(0.0, true)); // far from 7, flat region
        params.insert("sigma".into(), bounded(0.3, 1e-6, f64::INFINITY));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let data = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        let r = fit(&graph, vec![data], &default_options()).expect("fit runs");
        assert!(
            !r.success,
            "a collapsed-peak fit (R²<0, amplitude≈0) must be flagged unsuccessful, got success=true r2={}",
            r.r_squared
        );
        assert!(
            r.message.contains("degenerate_fit"),
            "expected degenerate_fit message, got {:?}",
            r.message
        );
    }

    // ── VarPro rejects dataset_index scoping (it is not scope-aware) ──────────

    #[test]
    fn varpro_rejects_dataset_index_scoped_graph() {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(1.0, true));
        params.insert("center".into(), make_param(0.0, true));
        params.insert("sigma".into(), make_param(1.0, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: Some(0), // per-dataset local node
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let x: Vec<f64> = (0..10).map(|i| i as f64).collect();
        let data = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y: vec![0.0; 10],
            sigma: None,
            label: None,
        };
        let opts = FitOptionsSpec {
            solver: "varpro".into(),
            ..default_options()
        };
        let err = fit(&graph, vec![data], &opts);
        assert!(
            err.is_err(),
            "varpro must reject dataset_index-scoped graphs rather than mis-scope them"
        );
        let msg = format!("{:?}", err.unwrap_err());
        assert!(
            msg.contains("dataset_index"),
            "error should mention dataset_index, got {msg}"
        );
    }

    // ── CX-VPE-01: VarPro routing must honour `Parameter.expr` ties ───────────

    /// Two unbounded, separable Gaussians. When `tie` is set, `g2.sigma` is tied
    /// to `g1.sigma` via `Parameter.expr` (no `expr_edge`), so the graph is
    /// VarPro-eligible on every axis EXCEPT the tie — isolating the tie check.
    fn two_gaussian_param_expr_graph(tie: bool) -> FitGraphSpec {
        let mut g1: HashMap<String, ParameterSpec> = HashMap::new();
        g1.insert("amplitude".into(), make_param(5.0, true));
        g1.insert("center".into(), make_param(-1.0, true));
        g1.insert("sigma".into(), make_param(0.5, true));

        let mut g2: HashMap<String, ParameterSpec> = HashMap::new();
        g2.insert("amplitude".into(), make_param(3.0, true));
        g2.insert("center".into(), make_param(1.5, true));
        // tied => vary=false and value derived from g1.sigma via Parameter.expr.
        let mut sig2 = make_param(0.5, !tie);
        if tie {
            sig2.expr = Some("g1.sigma".into());
        }
        g2.insert("sigma".into(), sig2);

        FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![
                ModelNodeSpec {
                    id: "g1".into(),
                    model_type: ModelTypeStr::Gaussian,
                    dataset_index: None,
                    parameters: g1,
                },
                ModelNodeSpec {
                    id: "g2".into(),
                    model_type: ModelTypeStr::Gaussian,
                    dataset_index: None,
                    parameters: g2,
                },
            ],
            expr_edges: vec![],
        }
    }

    #[test]
    fn graph_prefers_varpro_false_for_param_expr_tie() {
        // CX-VPE-01 regression: an otherwise VarPro-eligible (separable, single
        // dataset, all nonlinear params unbounded) graph whose ONLY tie lives in
        // `Parameter.expr` must NOT auto-route to VarPro (which would silently drop
        // the tie). Before the fix this returned true (only `expr_edges` checked).
        let graph = two_gaussian_param_expr_graph(true);
        assert!(
            !graph_prefers_varpro(&graph),
            "a Parameter.expr tie must disqualify the graph from VarPro auto-routing"
        );
    }

    #[test]
    fn graph_prefers_varpro_true_for_untied_unbounded_separable() {
        // Positive control: the SAME unbounded separable graph WITHOUT a tie stays
        // VarPro-eligible — the fix must not over-reject untied graphs.
        let graph = two_gaussian_param_expr_graph(false);
        assert!(
            graph_prefers_varpro(&graph),
            "an untied, unbounded, separable graph must remain VarPro-eligible"
        );
    }

    #[test]
    fn auto_and_explicit_varpro_agree_on_success() {
        // REGRESSION. The `Solver::Auto if graph_prefers_varpro` arm used to call
        // `solve_varpro_path` raw while `Solver::Varpro` wrapped it in
        // `finalize_varpro_result`. The guards live in that wrapper, so on a
        // VarPro-eligible graph the two solver strings could disagree about
        // `success` for identical data — with `auto`, the documented default,
        // being the one that reported the false success.
        //
        // The assertion is deliberately on AGREEMENT rather than on a specific
        // verdict: whether this particular fit converges is not the point, and
        // pinning it would make the test brittle to solver tuning. What must
        // hold is that routing does not change the answer.
        let graph = two_gaussian_param_expr_graph(false);
        assert!(
            graph_prefers_varpro(&graph),
            "fixture must actually route auto -> varpro, or this proves nothing"
        );

        let n = 64usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -3.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        // Data two Gaussians cannot represent: a high-frequency oscillation.
        // The fit converges but to a poor optimum, which is what puts r² under
        // OFF_DOMAIN_R2_FLOOR (0.5) so `apply_postfit_guards` actually engages.
        // A clean, well-fit dataset makes this test VACUOUS — the guard's own
        // r²-quality escape short-circuits it and both paths trivially agree.
        let y: Vec<f64> = x.iter().map(|xi| (20.0 * xi).sin()).collect();
        let data = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        let opts = |s: &str| FitOptionsSpec {
            solver: s.to_string(),
            ..default_options()
        };

        let r_auto = fit(&graph, vec![data.clone()], &opts("auto")).expect("auto fit");
        let r_varpro = fit(&graph, vec![data], &opts("varpro")).expect("varpro fit");

        assert_eq!(
            r_auto.success, r_varpro.success,
            "auto and explicit varpro must reach the same success verdict on the \
             same graph and data; a difference means one path skipped the \
             post-fit guards in `finalize_varpro_result`"
        );
        assert_relative_eq!(r_auto.chi2, r_varpro.chi2, max_relative = 1e-12);
    }

    #[test]
    fn varpro_explicit_rejects_param_expr_tie() {
        // Explicit solver="varpro" on a `Parameter.expr`-tied graph must return the
        // tied-params-unsupported error, NOT a silently-wrong success.
        let graph = two_gaussian_param_expr_graph(true);
        let n = 64usize;
        let x: Vec<f64> = (0..n)
            .map(|i| -3.0 + 6.0 * i as f64 / (n - 1) as f64)
            .collect();
        let y: Vec<f64> = x
            .iter()
            .map(|&xi| gaussian(xi, 5.0, -1.0, 0.5) + gaussian(xi, 3.0, 1.5, 0.5))
            .collect();
        let data = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        let opts = FitOptionsSpec {
            solver: "varpro".into(),
            ..default_options()
        };
        let msg = format!(
            "{}",
            fit(&graph, vec![data], &opts).expect_err(
                "solver='varpro' with a Parameter.expr tie must error, not drop the tie"
            )
        );
        assert!(
            msg.contains("tied parameters") || msg.contains("Parameter.expr"),
            "error should name the tied-params limitation, got: {msg}"
        );
    }

    /// S1 regression: a node with an UNBOUNDED amplitude but a BOUNDED nonlinear
    /// param (sigma) must NOT auto-route to VarPro (which ignores bounds). The old
    /// positional `.skip(1)` over a HashMap could skip sigma instead of amplitude
    /// and wrongly green-light VarPro — a silent bound violation.
    #[test]
    fn graph_with_unbounded_amplitude_but_bounded_sigma_is_not_varpro() {
        let bounded_sigma = ParameterSpec {
            value: 0.6,
            min: 0.1,
            max: 2.0,
            vary: true,
            expr: None,
            scale: None,
        };
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(4.0, true)); // unbounded linear coeff
        params.insert("center".into(), make_param(0.0, true)); // unbounded
        params.insert("sigma".into(), bounded_sigma); // BOUNDED nonlinear param
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        assert!(
            !graph_prefers_varpro(&graph),
            "a bounded nonlinear param must block VarPro auto-routing (bounds would be ignored)"
        );
    }

    // ------------------------------------------------------------------
    // A2 follow-up: typed SolverError variants reach the boundary
    // ------------------------------------------------------------------

    /// Build a minimal 1-D Gaussian graph + dataset for VarPro-rejection tests.
    fn make_varpro_inputs() -> (FitGraphSpec, MeasurementSpec) {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(1.0, true));
        params.insert("center".into(), make_param(0.0, true));
        params.insert("sigma".into(), make_param(1.0, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let x: Vec<f64> = (0..10).map(|i| i as f64 * 0.1).collect();
        let y: Vec<f64> = x.iter().map(|&xi| gaussian(xi, 1.0, 0.0, 1.0)).collect();
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        (graph, dataset)
    }

    /// VarPro with expr_edges must surface the typed
    /// `SolverError::VarproExprEdgesUnsupported` variant via the CoreError
    /// boundary conversion.
    #[test]
    fn varpro_with_expr_edges_emits_solver_error_variant() {
        use spectrafit_types::ExprEdge;

        let (mut graph, dataset) = make_varpro_inputs();
        // Tie one parameter to force the expr_edge check to fail.
        graph.nodes[0]
            .parameters
            .insert("amplitude".to_string(), make_param(1.0, false));
        graph.expr_edges.push(ExprEdge {
            target_node: "g1".to_string(),
            target_param: "amplitude".to_string(),
            expression: "2.0".to_string(),
        });

        let mut options = default_options();
        options.solver = "varpro".to_string();
        let err = fit(&graph, vec![dataset], &options).unwrap_err();
        let expected: CoreError = SolverError::VarproExprEdgesUnsupported.into();
        assert_eq!(format!("{err}"), format!("{expected}"));
    }

    /// VarPro with a `dataset_index`-scoped node must surface the typed
    /// `SolverError::VarproDatasetScopingUnsupported` variant.
    #[test]
    fn varpro_with_dataset_index_emits_solver_error_variant() {
        let (mut graph, dataset) = make_varpro_inputs();
        graph.nodes[0].dataset_index = Some(0);

        let mut options = default_options();
        options.solver = "varpro".to_string();
        let err = fit(&graph, vec![dataset], &options).unwrap_err();
        let expected: CoreError = SolverError::VarproDatasetScopingUnsupported.into();
        assert_eq!(format!("{err}"), format!("{expected}"));
    }

    // ------------------------------------------------------------------
    // Coverage-gap follow-up: `Solver::Varpro` on a non-separable graph
    // ------------------------------------------------------------------

    /// Explicit `solver="varpro"` on a graph that `spectrafit_varpro::is_separable`
    /// rejects (e.g. `true_voigt`, which is not in `SEPARABLE_MODEL_TYPES`) must
    /// surface `SolverError::VarproNotSeparable`, not silently fall back to
    /// another solver or panic.
    #[test]
    fn fit_rejects_varpro_on_non_separable_graph() {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(1.0, true));
        params.insert("center".into(), make_param(0.0, true));
        params.insert("sigma".into(), make_param(1.0, true));
        params.insert("gamma".into(), make_param(1.0, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "v1".into(),
                model_type: ModelTypeStr::TrueVoigt,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let x: Vec<f64> = (0..10).map(|i| i as f64 * 0.1).collect();
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y: vec![0.0; 10],
            sigma: None,
            label: None,
        };
        let opts = FitOptionsSpec {
            solver: "varpro".to_string(),
            ..default_options()
        };
        let err = fit(&graph, vec![dataset], &opts).expect_err(
            "solver='varpro' on a non-separable graph (true_voigt) must error, not fall back",
        );
        let expected: CoreError = SolverError::VarproNotSeparable.into();
        assert_eq!(format!("{err}"), format!("{expected}"));
    }

    // ------------------------------------------------------------------
    // Coverage-gap follow-up: `Solver::Global`'s max_iterations → DE-generation
    // budget mapping in `dispatch_solver`
    // ------------------------------------------------------------------

    /// A cheap single-parameter constant fit through `solver="global"`, used to
    /// exercise `dispatch_solver`'s DE-generation-budget mapping without the
    /// cost of a multi-peak DE search.
    fn constant_global_fit_inputs() -> (FitGraphSpec, MeasurementSpec) {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert(
            "c".into(),
            ParameterSpec {
                value: 0.0,
                min: -10.0,
                max: 10.0,
                vary: true,
                expr: None,
                scale: None,
            },
        );
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "const1".into(),
                model_type: ModelTypeStr::Constant,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        let n = 10usize;
        let x: Vec<f64> = (0..n).map(|i| i as f64).collect();
        let y = vec![3.0; n];
        let dataset = MeasurementSpec {
            schema_version: None,
            x: vec![x],
            y,
            sigma: None,
            label: None,
        };
        (graph, dataset)
    }

    /// `options.max_iterations > 0` must map onto `DeConfig::max_gen`, capping
    /// the DE search at that many generations — see the comment on
    /// `Solver::Global` in `dispatch_solver`.
    #[test]
    fn dispatch_global_solver_uses_max_iterations_as_de_generation_budget() {
        let (graph, dataset) = constant_global_fit_inputs();
        let opts = FitOptionsSpec {
            solver: "global".to_string(),
            max_iterations: 7,
            ..default_options()
        };
        let result = fit(&graph, vec![dataset], &opts).expect("global fit should not error");
        let gens = result
            .n_de_generations
            .expect("global solver must report n_de_generations");
        assert!(
            gens <= 7,
            "max_iterations=7 must cap the DE generation budget at 7, got {gens}"
        );
    }

    /// `options.max_iterations == 0` must fall back to `DeConfig::default().max_gen`
    /// (100) rather than passing a zero-generation budget through to DE.
    #[test]
    fn dispatch_global_solver_falls_back_to_default_generation_budget_when_max_iterations_is_zero()
    {
        let (graph, dataset) = constant_global_fit_inputs();
        let opts = FitOptionsSpec {
            solver: "global".to_string(),
            max_iterations: 0,
            ..default_options()
        };
        let result = fit(&graph, vec![dataset], &opts)
            .expect("global fit with max_iterations=0 must still run, using the DE default");
        let gens = result
            .n_de_generations
            .expect("global solver must report n_de_generations");
        assert!(
            gens <= DeConfig::default().max_gen as u64,
            "max_iterations=0 must fall back to DeConfig::default().max_gen (100), got {gens}"
        );
    }

    // ------------------------------------------------------------------
    // Coverage-gap follow-up: `run_lm_solve`'s tolerance/delta0/max_delta
    // fallback and override branches
    // ------------------------------------------------------------------

    /// `options.tolerance <= 0.0` must fall back to the built-in `1e-8` default
    /// for the shared `tol` used by the faer/trust-region drivers (the `else`
    /// arm of the top-level `tol` binding), AND must skip `LevenbergMarquardt::with_tol`
    /// for the `lm-legacy` oracle (its own, separate `options.tolerance > 0.0`
    /// check) rather than calling it with a non-positive tolerance.
    #[test]
    fn run_lm_solve_falls_back_to_default_tolerance_when_non_positive() {
        let (graph, dataset) = gaussian_fit_inputs(None);
        let opts = FitOptionsSpec {
            solver: "lm-legacy".to_string(),
            tolerance: 0.0,
            ..default_options()
        };
        let result = fit(&graph, vec![dataset], &opts)
            .expect("lm-legacy fit with tolerance<=0.0 should still run on the 1e-8 default");
        assert!(
            result.success,
            "lm-legacy should converge on a clean Gaussian even via the default tolerance; \
             message: {}",
            result.message
        );
    }

    /// An explicit `delta0`/`max_delta` must override the trust-region driver's
    /// own defaults (the `if let Some(..)` bodies in `run_lm_solve`), not be
    /// silently ignored.
    #[test]
    fn run_lm_solve_applies_explicit_delta0_and_max_delta_overrides() {
        let (graph, dataset) = gaussian_fit_inputs(None);
        let opts = FitOptionsSpec {
            solver: "dogleg".to_string(),
            delta0: Some(0.5),
            max_delta: Some(50.0),
            ..default_options()
        };
        let result = fit(&graph, vec![dataset], &opts)
            .expect("dogleg fit with explicit delta0/max_delta should not error");
        assert!(
            result.success,
            "dogleg should converge on a clean Gaussian with an explicit \
             delta0/max_delta override; message: {}",
            result.message
        );
    }

    // ------------------------------------------------------------------
    // Coverage-gap follow-up: `collect_free_param_specs`'s malformed-lookup
    // error path
    // ------------------------------------------------------------------

    /// A `free_keys` entry naming a parameter the node itself does not carry
    /// (e.g. a stale/malformed key) must surface a `SolverError::Dispatch`
    /// naming the missing parameter and node, not panic on an unwrap.
    #[test]
    fn collect_free_param_specs_errors_when_free_key_param_missing_from_node() {
        let mut params: HashMap<String, ParameterSpec> = HashMap::new();
        params.insert("amplitude".into(), make_param(1.0, true));
        let graph = FitGraphSpec {
            schema_version: "0.1".into(),
            nodes: vec![ModelNodeSpec {
                id: "g1".into(),
                model_type: ModelTypeStr::Gaussian,
                dataset_index: None,
                parameters: params,
            }],
            expr_edges: vec![],
        };
        // "sigma" is a plausible free key for a Gaussian node, but this node's
        // own parameter map only has "amplitude" — collect_free_param_specs must
        // report the mismatch rather than panicking on a missing-key unwrap.
        let free_keys = vec!["g1.sigma".to_string()];
        let err = collect_free_param_specs(&graph, &free_keys)
            .expect_err("a free key naming a param absent from its node must error");
        let msg = format!("{err}");
        assert!(
            msg.contains("sigma") && msg.contains("g1"),
            "error should name the missing param and its node, got: {msg}"
        );
    }

    // ------------------------------------------------------------------
    // Coverage-gap follow-up: termination-reason mapping helpers
    // ------------------------------------------------------------------

    /// `faer_termination_str` must map every `spectrafit_levenberg_marquardt::Termination`
    /// variant to its stable snake_case string — including `NumericalError`,
    /// which no end-to-end fit in this suite naturally triggers.
    #[test]
    fn faer_termination_str_maps_numerical_error() {
        assert_eq!(
            faer_termination_str(spectrafit_levenberg_marquardt::Termination::NumericalError),
            "numerical_error"
        );
    }

    /// `map_termination` must translate every `levenberg_marquardt::TerminationReason`
    /// variant (the `lm-legacy` oracle's own type) onto our stable
    /// `spectrafit_types::TerminationReason`, independent of any particular fit
    /// reaching that termination naturally.
    #[test]
    fn map_termination_covers_every_lm_legacy_variant() {
        use levenberg_marquardt::TerminationReason as L;

        let cases: Vec<(L, TerminationReason)> = vec![
            (L::ResidualsZero, TerminationReason::ResidualsZero),
            (L::Orthogonal, TerminationReason::Orthogonal),
            (L::LostPatience, TerminationReason::MaxIterations),
            (
                L::NoImprovementPossible("test"),
                TerminationReason::NoImprovementPossible,
            ),
            (L::NoParameters, TerminationReason::NoParameters),
            (L::NoResiduals, TerminationReason::NoResiduals),
            (
                L::WrongDimensions("test"),
                TerminationReason::WrongDimensions,
            ),
            (L::Numerical("test"), TerminationReason::NumericalError),
            (L::User("test"), TerminationReason::UserCancelled),
        ];
        for (input, expected) in cases {
            let mapped = map_termination(&input);
            assert_eq!(
                mapped, expected,
                "map_termination({input:?}) should map to {expected:?}, got {mapped:?}"
            );
        }
    }
}