nereids_pipeline/

pipeline.rs

//! Single-spectrum analysis pipeline.
//!
//! Orchestrates the full analysis chain for a single transmission spectrum:
//! ENDF loading → cross-section calculation → broadening → fitting.
//!
//! This is the building block for the spatial mapping pipeline.
//!
//! Uses `UnifiedFitConfig` with `SolverConfig` and typed `InputData` variants.

use std::fmt;
use std::sync::Arc;

use nereids_core::constants::{EV_TO_JOULES, NEUTRON_MASS_KG};
use nereids_endf::resonance::ResonanceData;
use nereids_fitting::joint_poisson::{self, JointPoissonFitConfig, JointPoissonObjective};
use nereids_fitting::lm::{self, FitModel, LmConfig, LmResult};
use nereids_fitting::parameters::{FitParameter, ParameterSet};
use nereids_fitting::poisson::{self, PoissonConfig};
use nereids_fitting::transmission_model::{
    EnergyScaleTransmissionModel, NormalizedTransmissionModel, PrecomputedTransmissionModel,
    TransmissionFitModel,
};
use nereids_physics::resolution::ResolutionFunction;
use nereids_physics::transmission::InstrumentParams;

use crate::error::PipelineError;

/// Working-grid σ + its layout (issue #608): Doppler-broadened σ on the working
/// grid paired with the data-index map back to the data grid.  Injected by
/// [`spatial_map_typed`] into [`UnifiedFitConfig`] for the precomputed
/// Gaussian-resolution path.
type PrecomputedWorkXs = (
    Arc<Vec<Vec<f64>>>,
    Arc<nereids_physics::transmission::WorkingGridLayout>,
);

/// SAMMY-style normalization and background configuration.
///
/// When enabled, the transmission model becomes:
///   T_out(E) = Anorm × T_inner(E) + BackA + BackB / √E + BackC × √E
///            + BackD × exp(−BackF / √E)
///
/// The first 4 background parameters (Anorm, BackA, BackB, BackC) are always
/// available.  The exponential tail (BackD, BackF) is optional and disabled
/// by default (`fit_back_d = false`, `fit_back_f = false`).
///
/// ## SAMMY Reference
/// SAMMY manual Sec III.E.2 — NORMAlization and BACKGround cards.
/// SAMMY fits up to 6 background terms; we implement all 6.
#[derive(Debug, Clone)]
pub struct BackgroundConfig {
    /// Initial value for the normalization factor (default 1.0).
    pub anorm_init: f64,
    /// Initial value for the constant background (default 0.0).
    pub back_a_init: f64,
    /// Initial value for the 1/√E background term (default 0.0).
    pub back_b_init: f64,
    /// Initial value for the √E background term (default 0.0).
    pub back_c_init: f64,
    /// Initial value for the exponential amplitude (default 0.01).
    ///
    /// Must be > 0 when `fit_back_f` is true, otherwise the Jacobian
    /// column for BackF is identically zero and BackF cannot be learned.
    pub back_d_init: f64,
    /// Initial value for the exponential decay constant (default 1.0).
    ///
    /// Units: √eV.  Must be > 0 when `fit_back_d` is true, otherwise
    /// BackD is indistinguishable from BackA (both become constants).
    pub back_f_init: f64,
    /// Whether Anorm is free (true) or fixed (false).
    pub fit_anorm: bool,
    /// Whether BackA is free (true) or fixed (false).
    pub fit_back_a: bool,
    /// Whether BackB is free (true) or fixed (false).
    pub fit_back_b: bool,
    /// Whether BackC is free (true) or fixed (false).
    pub fit_back_c: bool,
    /// Whether BackD (exponential amplitude) is free (true) or fixed (false).
    pub fit_back_d: bool,
    /// Whether BackF (exponential decay constant) is free (true) or fixed (false).
    pub fit_back_f: bool,
}

impl Default for BackgroundConfig {
    fn default() -> Self {
        Self {
            anorm_init: 1.0,
            back_a_init: 0.0,
            back_b_init: 0.0,
            back_c_init: 0.0,
            back_d_init: 0.01,
            back_f_init: 1.0,
            fit_anorm: true,
            fit_back_a: true,
            fit_back_b: true,
            fit_back_c: true,
            fit_back_d: false,
            fit_back_f: false,
        }
    }
}

/// Indices of SAMMY background parameters in the full parameter vector.
///
/// Replaces the previous 4-tuple representation and supports the optional
/// exponential tail terms (BackD, BackF).
#[derive(Debug, Clone, Copy)]
struct BackgroundIndices {
    anorm: usize,
    back_a: usize,
    back_b: usize,
    back_c: usize,
    back_d: Option<usize>,
    back_f: Option<usize>,
}

// ── New typed pipeline API (Phase 0) ─────────────────────────────────────

/// Typed input data — makes the data format explicit at the API boundary.
///
/// The two variants carry genuinely different data:
/// - **Counts**: raw detector counts + open beam → Poisson statistics native
/// - **Transmission**: normalized T = sample/open_beam + uncertainty → Gaussian statistics
///
/// `spatial_map_typed` dispatches to the correct fitting engine based on
/// which variant is provided.  This eliminates the overloaded positional
/// arguments that caused silent misinterpretation in the old API.
#[derive(Debug, Clone)]
pub enum InputData {
    /// Pre-normalized transmission with Gaussian uncertainties.
    ///
    /// Use with LM (default) or Poisson KL (opt-in for low-count T data).
    Transmission {
        /// Measured transmission T(E), values typically in [0, 2].
        transmission: Vec<f64>,
        /// Per-bin uncertainty σ_T(E).
        uncertainty: Vec<f64>,
    },
    /// Raw detector counts with open beam reference.
    ///
    /// Always uses Poisson KL (statistically optimal for count data).
    /// The fitting engine works directly on counts — no information-losing
    /// normalization to transmission.
    Counts {
        /// Sample counts per energy bin.
        sample_counts: Vec<f64>,
        /// Open-beam counts per energy bin (normalization reference).
        open_beam_counts: Vec<f64>,
    },
    /// Counts with pre-estimated nuisance parameters (power users).
    ///
    /// Use when you want to inspect or override the flux estimate before fitting.
    CountsWithNuisance {
        /// Sample counts per energy bin.
        sample_counts: Vec<f64>,
        /// Estimated flux spectrum (from open beam spatial average).
        flux: Vec<f64>,
        /// Estimated detector background spectrum.
        background: Vec<f64>,
    },
}

impl InputData {
    /// Number of energy bins.
    pub fn n_energies(&self) -> usize {
        match self {
            Self::Transmission { transmission, .. } => transmission.len(),
            Self::Counts { sample_counts, .. } => sample_counts.len(),
            Self::CountsWithNuisance { sample_counts, .. } => sample_counts.len(),
        }
    }

    /// Whether this is count data (Poisson-native).
    pub fn is_counts(&self) -> bool {
        matches!(self, Self::Counts { .. } | Self::CountsWithNuisance { .. })
    }
}

/// Solver-specific configuration.
///
/// Carries the full solver config inside each variant, making invalid
/// combinations unrepresentable.
#[derive(Debug, Clone, Default)]
pub enum SolverConfig {
    /// Levenberg-Marquardt chi-squared minimizer (transmission path).
    LevenbergMarquardt(LmConfig),
    /// Poisson-KL counts-domain fitter.
    ///
    /// For **counts** inputs this dispatches to the joint-Poisson profile-
    /// binomial-deviance path (`joint_poisson_fit`) validated in memo 35
    /// §P1/§P2 and memo 38.  Uses an explicit `c = Q_s/Q_ob` from
    /// `CountsBackgroundConfig::c` and reports `D/(n − k)` as the primary
    /// GOF.  Stage-1 damped Fisher + optional Nelder-Mead polish (see
    /// [`nereids_fitting::joint_poisson::JointPoissonFitConfig`]).
    ///
    /// For **transmission** inputs this dispatches to Poisson NLL on the
    /// transmission values directly (legacy path, unchanged).  The payload
    /// `PoissonConfig` carries `max_iter` common to both dispatches.
    PoissonKL(PoissonConfig),
    /// Automatic: Counts → PoissonKL (counts-domain joint-Poisson),
    /// Transmission → LM.
    #[default]
    Auto,
}

/// Background model for the counts fitting engine.
///
/// In the counts domain, the forward model is:
///   Y(E) = α₁ · [Φ(E) · exp(-Σ nᵢσᵢ(E))] + α₂ · B(E)
///
/// where Φ(E) is the incident flux and B(E) is detector/gamma background.
/// The reference Φ(E) / B(E) spectra are supplied by the caller or by
/// spatial pre-processing; this config only controls the fitted scale factors.
///
/// Important distinction:
/// - This is a detector-space counts background model `B(E)`.
/// - It is NOT the same as the transmission-lift background used by
///   `BackgroundConfig`, which models additive uplift of the apparent
///   transmission curve (for example gamma-tail structure that pushes
///   transmission upward).
///
/// For VENUS MCP/TPX event detectors, the current working assumption is:
/// - raw/open-beam is the correct normalization baseline
/// - dark-current / CCD-style electronic offset is not modeled
/// - rare ghost counts may exist at the hardware level, but are currently
///   treated as negligible unless a detector-background reference spectrum
///   is explicitly provided
///
/// This is structurally different from the transmission background model
/// ([`BackgroundConfig`]) because:
/// - Φ and B are reference spectra, not fitted per pixel
/// - α₁ and α₂ are optional per-pixel scale corrections
/// - All terms are non-negative (required for valid Poisson NLL)
#[derive(Debug, Clone)]
pub struct CountsBackgroundConfig {
    /// **Research-only.** Initial α₁ flux-scale value used by
    /// [`crate::pipeline::evaluate_jacobian_and_fisher`] (Fisher-info
    /// research helper, Epic #394).  Not honoured by the production fit
    /// path; `SolverConfig::PoissonKL` profiles the flux via `λ̂` and
    /// rejects alpha fitting.
    pub alpha_1_init: f64,
    /// **Research-only.** Initial α₂ detector-bg-scale value, same
    /// provisions as `alpha_1_init`.
    pub alpha_2_init: f64,
    /// **Research-only.** Fit α₁ flag — only consumed by
    /// `evaluate_jacobian_and_fisher`.  Passing `true` through
    /// `SolverConfig::PoissonKL` on counts input yields an error.
    pub fit_alpha_1: bool,
    /// **Research-only.** Fit α₂ flag — see `fit_alpha_1`.
    pub fit_alpha_2: bool,
    /// Proton-charge ratio `c = Q_s / Q_ob` for the counts-KL solver
    /// (memo 35 §P1.3 — "make `c` a first-class API parameter").
    ///
    /// Default `1.0`, correct only when the caller has already PC-
    /// normalized the open-beam counts so that `flux = c · O`.  For the
    /// counts-KL dispatch (`SolverConfig::PoissonKL` on counts input),
    /// set this to the actual `Q_sample / Q_open_beam` ratio and pass
    /// raw open-beam counts — the joint-Poisson solver will profile
    /// `λ̂` itself.  Ignored by LM paths.
    pub c: f64,
}

impl Default for CountsBackgroundConfig {
    fn default() -> Self {
        Self {
            alpha_1_init: 1.0,
            alpha_2_init: 1.0,
            fit_alpha_1: false,
            fit_alpha_2: false,
            c: 1.0,
        }
    }
}

// ── Phase 2: UnifiedFitConfig + fit_spectrum_typed ───────────────────────

/// Unified fit configuration for all data types and solvers.
///
/// Carries both transmission and counts background configs, and uses
/// [`SolverConfig`] (which embeds solver-specific tuning).
#[derive(Debug, Clone)]
pub struct UnifiedFitConfig {
    // ── Physics (shared by both engines) ──
    energies: Vec<f64>,
    resonance_data: Vec<ResonanceData>,
    isotope_names: Vec<String>,
    temperature_k: f64,
    resolution: Option<ResolutionFunction>,
    initial_densities: Vec<f64>,
    fit_temperature: bool,
    compute_covariance: bool,

    // ── Solver ──
    solver: SolverConfig,

    // ── Background models (engine-specific) ──
    /// SAMMY-style background for the transmission engine.
    transmission_background: Option<BackgroundConfig>,
    /// Counts-domain background for the counts engine.
    counts_background: Option<CountsBackgroundConfig>,

    // ── Joint-Poisson solver knobs (counts path only) ──
    /// When `Some(false)`, the counts-KL dispatch disables Nelder-Mead polish
    /// (stage-1 damped Fisher only).  When `Some(true)`, polish is forced on
    /// regardless of context.  When `None`, the dispatcher picks a default:
    /// polish on for single-spectrum fits, off for per-pixel spatial maps
    /// (memo 38 §6 recommendation — 17 min polish per pixel is untenable).
    counts_enable_polish: Option<bool>,

    // ── Precomputed caches (injected by spatial_map_typed) ──
    precomputed_cross_sections: Option<Arc<Vec<Vec<f64>>>>,
    /// Doppler-broadened σ on the **working grid** + the working-grid layout,
    /// injected by [`spatial_map_typed`] for the fixed-calibration /
    /// fixed-temperature precomputed path (issue #608).
    ///
    /// When a Gaussian resolution function is active, the working grid is the
    /// auxiliary extended grid (boundary extension + resonance fine-structure);
    /// storing σ there lets each per-pixel [`PrecomputedTransmissionModel`]
    /// apply Beer-Lambert + resolution on the working grid and extract the data
    /// points last — matching `forward_model`.  For tabulated / no resolution
    /// the working grid IS the data grid and this is `None` (the model uses the
    /// data-grid `precomputed_cross_sections` directly, preserving the cubature
    /// / scalar surrogate fast paths).
    ///
    /// `precomputed_cross_sections` still carries the **data-grid** σ for the
    /// surrogate-plan builders and shape validation; this field is the separate
    /// working-grid copy consumed only by `build_transmission_model`'s
    /// precomputed branch.
    precomputed_work_cross_sections: Option<PrecomputedWorkXs>,
    precomputed_base_xs: Option<Arc<Vec<Vec<f64>>>>,
    /// Resolution broadening plan built once for `(energies, resolution)`.
    ///
    /// Populated by [`spatial_map_typed`] when the data grid is shared
    /// across every pixel, so each per-pixel fit reuses a single
    /// hoisted TOF / kernel-interpolation / bracket table instead of
    /// rebuilding it every broadening call.  `None` ⇒ the fit-model
    /// layer falls back to the per-call broadening path with byte-
    /// identical output.
    precomputed_resolution_plan: Option<Arc<nereids_physics::resolution::ResolutionPlan>>,
    /// Sparse empirical cubature plan built once for `(energies,
    /// isotope_set, density_box)` — when present and the dispatch
    /// conditions hold (k ≥ 2, fixed calibration, fixed temperature),
    /// the fit model replaces `exp(-Σ n σ) + apply_resolution` with
    /// `cubature.forward_and_jacobian(n)` directly.  See epic #472.
    ///
    /// `None` ⇒ existing path.  The cubature is advisory — if its
    /// target grid doesn't match `energies` or if `k == 1` or
    /// temperature/energy-scale fitting is active, the fit model
    /// silently falls back to the exact path.
    precomputed_sparse_cubature_plan:
        Option<Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>>,
    /// Scalar (k = 1) surrogate plan for the grouped-Hf / single-
    /// isotope path.  Parallel to `precomputed_sparse_cubature_plan`
    /// but dispatches at `k == 1`.  Epic #472, PR #475.
    precomputed_sparse_scalar_plan: Option<Arc<nereids_physics::surrogate::ScalarSurrogatePlan>>,

    // ── Energy-scale calibration (SAMMY TZERO equivalent) ──
    /// When true, fit t₀ (μs) and L_scale (dimensionless) parameters.
    /// These adjust the energy axis during fitting:
    ///   E_corr = (TOF_FACTOR * L * L_scale / (t_nom - t₀))²
    fit_energy_scale: bool,
    /// Initial t₀ value in microseconds (default 0.0).
    t0_init_us: f64,
    /// Initial L_scale value (dimensionless, default 1.0).
    l_scale_init: f64,
    /// Flight path in meters for TOF↔energy conversion (default from resolution or 25.0).
    flight_path_m: f64,
    /// Method for the t0 / L_scale Jacobian columns. `None` defers to
    /// `EnergyScaleJacobianMethod::from_env`: it uses the
    /// `NEREIDS_TZERO_JACOBIAN` env var when set, and otherwise falls
    /// back to `PartialGal` (the default since issue #489).
    /// `Some(_)` bypasses both the env var and the default.
    tzero_jacobian_method: Option<nereids_fitting::transmission_model::EnergyScaleJacobianMethod>,

    // ── Fit energy range restriction (SAMMY EMIN/EMAX equivalent) ──
    /// User-specified fit-energy-range restriction.  When `Some((min,
    /// max))`, the configured `energies` grid is expected to extend
    /// beyond `[min, max]` by a kernel-margin (~5×FWHM); the GUI / pre-
    /// processing layer slices the input data to that extended grid
    /// before constructing this config.  The solver cost functions
    /// (LM transmission and joint-Poisson PBD) mask residuals outside
    /// `[min, max]` to zero so resonance contributions just inside the
    /// boundaries are correctly broadened.
    ///
    /// `None` (default) = full grid, no masking.  Reduced χ² / dof
    /// counts are computed against the active bin count when masking
    /// is in effect.  SAMMY equivalent: the `EMIN` / `EMAX` analysis
    /// limits (INPut-file card set 2, manual Table VI A.1); the kernel
    /// margin follows the same endpoint-extension principle as SAMMY's
    /// auxiliary grid (general construction Sec. III.A.2(c); the
    /// quantitative `[Emin − Wmin, Emax + Wmax]` statement is in the
    /// Leal-Hwang procedure, Sec. III.B.2), with a deliberately
    /// conservative 5×FWHM margin.
    fit_energy_range: Option<(f64, f64)>,

    // ── Isotope group mapping (optional) ──
    /// Maps member isotope index → density parameter index.
    /// `None` = identity mapping (one param per isotope, backward compat).
    pub(crate) density_indices: Option<Vec<usize>>,
    /// Fractional ratio per member isotope.
    /// `None` = all 1.0 (backward compat).
    pub(crate) density_ratios: Option<Vec<f64>>,
    /// Number of density parameters (groups or isotopes).
    /// `None` = `resonance_data.len()` (backward compat).
    n_density_params: Option<usize>,
}

impl UnifiedFitConfig {
    /// Construct a new config with validation.
    pub fn new(
        energies: Vec<f64>,
        resonance_data: Vec<ResonanceData>,
        isotope_names: Vec<String>,
        temperature_k: f64,
        resolution: Option<ResolutionFunction>,
        initial_densities: Vec<f64>,
    ) -> Result<Self, FitConfigError> {
        if energies.is_empty() {
            return Err(FitConfigError::EmptyEnergies);
        }
        if resonance_data.is_empty() {
            return Err(FitConfigError::EmptyResonanceData);
        }
        if initial_densities.len() != resonance_data.len() {
            return Err(FitConfigError::DensityCountMismatch {
                densities: initial_densities.len(),
                isotopes: resonance_data.len(),
            });
        }
        if isotope_names.len() != resonance_data.len() {
            return Err(FitConfigError::NameCountMismatch {
                names: isotope_names.len(),
                isotopes: resonance_data.len(),
            });
        }
        if !temperature_k.is_finite() {
            return Err(FitConfigError::NonFiniteTemperature(temperature_k));
        }
        if temperature_k < 0.0 {
            return Err(FitConfigError::NegativeTemperature(temperature_k));
        }
        Ok(Self {
            energies,
            resonance_data,
            isotope_names,
            temperature_k,
            resolution,
            initial_densities,
            fit_temperature: false,
            compute_covariance: true,
            solver: SolverConfig::Auto,
            transmission_background: None,
            counts_background: None,
            counts_enable_polish: None,
            precomputed_cross_sections: None,
            precomputed_work_cross_sections: None,
            precomputed_base_xs: None,
            precomputed_resolution_plan: None,
            precomputed_sparse_cubature_plan: None,
            precomputed_sparse_scalar_plan: None,
            fit_energy_scale: false,
            t0_init_us: 0.0,
            l_scale_init: 1.0,
            flight_path_m: 25.0,
            density_indices: None,
            density_ratios: None,
            n_density_params: None,
            tzero_jacobian_method: None,
            fit_energy_range: None,
        })
    }

    /// Override the method used for the t0 / L_scale Jacobian columns
    /// in `EnergyScaleTransmissionModel`. `None` (default) defers to
    /// the model's own default selection via
    /// `EnergyScaleJacobianMethod::from_env`: `PartialGal` since issue
    /// #489 unless overridden by `NEREIDS_TZERO_JACOBIAN`.
    #[must_use]
    pub fn with_tzero_jacobian_method(
        mut self,
        method: Option<nereids_fitting::transmission_model::EnergyScaleJacobianMethod>,
    ) -> Self {
        self.tzero_jacobian_method = method;
        self
    }

    // ── Builder methods ──

    #[must_use]
    pub fn with_solver(mut self, solver: SolverConfig) -> Self {
        self.solver = solver;
        self
    }

    #[must_use]
    pub fn with_fit_temperature(mut self, v: bool) -> Self {
        self.fit_temperature = v;
        self
    }

    #[must_use]
    pub fn with_compute_covariance(mut self, v: bool) -> Self {
        self.compute_covariance = v;
        self
    }

    /// Enable energy-scale fitting (SAMMY TZERO equivalent).
    ///
    /// Adds t₀ (μs) and L_scale (dimensionless) as fit parameters.
    /// These adjust the energy axis during fitting to correct for
    /// flight-path and timing-offset uncertainties.
    #[must_use]
    pub fn with_energy_scale(
        mut self,
        t0_init_us: f64,
        l_scale_init: f64,
        flight_path_m: f64,
    ) -> Self {
        self.fit_energy_scale = true;
        self.t0_init_us = t0_init_us;
        self.l_scale_init = l_scale_init;
        self.flight_path_m = flight_path_m;
        self
    }

    /// Restrict the fit cost function to bins inside `[min_eV, max_eV]`
    /// (SAMMY EMIN/EMAX equivalent).  The configured `energies` grid is
    /// expected to extend by ~5×FWHM beyond `[min, max]` on each side
    /// (the GUI / pre-processing layer handles this); the LM and
    /// joint-Poisson cost paths mask residuals outside `[min, max]` to
    /// zero so resonance broadening at the boundaries is correct.
    /// `None` (default) = full grid, no masking.
    ///
    /// Validates the range up-front so non-finite or reversed bounds
    /// (which would silently produce an empty active-bin mask and a
    /// deceptive fit) are rejected at config-build time rather than
    /// surfacing as a downstream solver error.
    pub fn with_fit_energy_range(
        mut self,
        range: Option<(f64, f64)>,
    ) -> Result<Self, FitConfigError> {
        if let Some((lo, hi)) = range {
            if !lo.is_finite() || !hi.is_finite() {
                return Err(FitConfigError::InvalidFitEnergyRange(
                    "fit_energy_range bounds must be finite",
                ));
            }
            if lo >= hi {
                return Err(FitConfigError::InvalidFitEnergyRange(
                    "fit_energy_range min must be strictly less than max",
                ));
            }
        }
        self.fit_energy_range = range;
        Ok(self)
    }

    #[must_use]
    pub fn with_transmission_background(mut self, bg: BackgroundConfig) -> Self {
        self.transmission_background = Some(bg);
        self
    }

    #[must_use]
    pub fn with_counts_background(mut self, bg: CountsBackgroundConfig) -> Self {
        self.counts_background = Some(bg);
        self
    }

    /// Override the Nelder-Mead polish flag for the counts-KL dispatch.
    /// `Some(true)` forces polish on, `Some(false)` forces it off, `None`
    /// (the default) lets the dispatcher pick (polish on for single-spectrum,
    /// off for spatial maps per memo 38 §6).
    #[must_use]
    pub fn with_counts_enable_polish(mut self, v: Option<bool>) -> Self {
        self.counts_enable_polish = v;
        self
    }

    #[must_use]
    pub fn with_precomputed_cross_sections(mut self, xs: Arc<Vec<Vec<f64>>>) -> Self {
        self.precomputed_cross_sections = Some(xs);
        // A new XS cache invalidates any prebuilt cubature — the
        // cubature's atoms encode the OLD σ stack, so reusing the
        // plan would silently produce wrong forward / Jacobian
        // values for the new σ (Codex round-3 P3 on PR #480).
        self.precomputed_sparse_cubature_plan = None;
        self.precomputed_sparse_scalar_plan = None;
        // A new data-grid σ also invalidates the working-grid σ copy
        // (issue #608): it was Doppler-broadened from the OLD σ on the
        // OLD grid layout, so reusing it would mix grids.
        self.precomputed_work_cross_sections = None;
        self
    }

    /// Attach the **working-grid** Doppler-broadened σ + its layout for the
    /// fixed-calibration / fixed-temperature precomputed path (issue #608).
    ///
    /// When set, `build_transmission_model` builds a
    /// [`PrecomputedTransmissionModel`] whose σ live on the working grid and
    /// whose `evaluate` / `analytical_jacobian` apply resolution on the working
    /// grid and extract the data points last — matching `forward_model`.  The
    /// data-grid `precomputed_cross_sections` must still be set (for the
    /// surrogate-plan builders and shape validation); for tabulated / no
    /// resolution the working grid equals the data grid and this is left
    /// `None`.
    ///
    /// The σ + layout are not validated here (this is an infallible builder
    /// setter); shape/consistency are checked once up front in
    /// `validate_precomputed_cross_sections`, which every public entry point
    /// (`fit_spectrum_typed`, `fit_transmission_poisson`, `spatial_map_typed`)
    /// calls before any forward-model build or per-pixel loop — mirroring how
    /// [`Self::with_precomputed_cross_sections`] is validated.
    #[must_use]
    pub fn with_precomputed_work_cross_sections(
        mut self,
        xs: Arc<Vec<Vec<f64>>>,
        layout: Arc<nereids_physics::transmission::WorkingGridLayout>,
    ) -> Self {
        self.precomputed_work_cross_sections = Some((xs, layout));
        self
    }

    #[must_use]
    pub fn with_precomputed_base_xs(mut self, xs: Arc<Vec<Vec<f64>>>) -> Self {
        self.precomputed_base_xs = Some(xs);
        // Base XS swap implies σ will be re-Doppler-broadened, so
        // the cubature's σ stack becomes stale.  Invalidate for
        // the same reason as `with_precomputed_cross_sections`.
        self.precomputed_sparse_cubature_plan = None;
        self.precomputed_sparse_scalar_plan = None;
        self
    }

    /// Attach a prebuilt resolution plan for the config's energy grid.
    ///
    /// The caller (typically `spatial_map_typed`) must ensure that
    /// `plan.target_energies()` equals `self.energies()`, otherwise
    /// the fit-model layer will return either a length-mismatch
    /// error or `ResolutionError::PlanGridMismatch` (for a different
    /// same-length grid) on the first broadening call.
    #[must_use]
    pub fn with_precomputed_resolution_plan(
        mut self,
        plan: Arc<nereids_physics::resolution::ResolutionPlan>,
    ) -> Self {
        self.precomputed_resolution_plan = Some(plan);
        self
    }

    /// Attach a prebuilt sparse empirical cubature plan for the
    /// config's energy grid + isotope set (see epic #472).
    ///
    /// The plan is advisory — the fit model falls back to the exact
    /// `ResolutionPlan` path when any of these guards fire:
    ///
    /// * `plan.target_energies() != self.energies()` (grid mismatch).
    /// * `plan.k() != n_density_params` (isotope-set mismatch).
    /// * `self.fit_temperature == true` (σ changes → atoms stale).
    /// * `self.fit_energy_scale == true` (grid changes → plan stale).
    /// * `n_density_params == 1` (scalar fast-path belongs to PR #475;
    ///   for now this path still falls back).
    ///
    /// Callers (typically `spatial_map_typed`) are responsible for
    /// ensuring the plan was built against compatible `sigmas` /
    /// `training_densities` / `jacobian_anchor`; the fit model cannot
    /// re-check those at dispatch time.
    #[must_use]
    pub fn with_precomputed_sparse_cubature_plan(
        mut self,
        plan: Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>,
    ) -> Self {
        self.precomputed_sparse_cubature_plan = Some(plan);
        self
    }

    /// Attach a prebuilt scalar (k = 1) surrogate plan.  Epic #472,
    /// PR #475.  Same invalidation discipline as the cubature: the
    /// plan is cleared on `with_groups` / `with_precomputed_cross_sections`
    /// / `with_precomputed_base_xs`, so a stale σ cannot silently
    /// dispatch.
    #[must_use]
    pub fn with_precomputed_sparse_scalar_plan(
        mut self,
        plan: Arc<nereids_physics::surrogate::ScalarSurrogatePlan>,
    ) -> Self {
        self.precomputed_sparse_scalar_plan = Some(plan);
        self
    }

    /// Configure isotope groups with ratio constraints.
    ///
    /// Each group binds multiple isotopes to one fitted density parameter.
    /// `groups` is a slice of `(IsotopeGroup, member_resonance_data)` pairs.
    /// `initial_densities` must have one entry per group.
    ///
    /// Replaces the existing per-isotope configuration with the expanded
    /// group mapping (flattened resonance_data + density_indices + density_ratios).
    pub fn with_groups(
        mut self,
        groups: &[(&nereids_core::types::IsotopeGroup, &[ResonanceData])],
        initial_densities: Vec<f64>,
    ) -> Result<Self, FitConfigError> {
        if groups.is_empty() {
            return Err(FitConfigError::EmptyResonanceData);
        }
        if initial_densities.len() != groups.len() {
            return Err(FitConfigError::DensityCountMismatch {
                densities: initial_densities.len(),
                isotopes: groups.len(),
            });
        }
        let mut all_resonance_data = Vec::new();
        let mut all_indices = Vec::new();
        let mut all_ratios = Vec::new();
        let mut names = Vec::new();
        for (g_idx, (group, rd_list)) in groups.iter().enumerate() {
            if rd_list.len() != group.n_members() {
                return Err(FitConfigError::GroupMemberCountMismatch {
                    group_name: group.name().to_string(),
                    rd_count: rd_list.len(),
                    member_count: group.n_members(),
                });
            }
            names.push(group.name().to_string());
            for ((isotope, ratio), rd) in group.members().iter().zip(rd_list.iter()) {
                // Validate that the ResonanceData matches the expected member isotope.
                if rd.isotope != *isotope {
                    return Err(FitConfigError::GroupMemberIsotopeMismatch {
                        group_name: group.name().to_string(),
                        expected_z: isotope.z(),
                        expected_a: isotope.a(),
                        got_z: rd.isotope.z(),
                        got_a: rd.isotope.a(),
                    });
                }
                all_resonance_data.push(rd.clone());
                all_indices.push(g_idx);
                all_ratios.push(*ratio);
            }
        }
        self.resonance_data = all_resonance_data;
        self.isotope_names = names;
        self.initial_densities = initial_densities;
        self.n_density_params = Some(groups.len());
        self.density_indices = Some(all_indices);
        self.density_ratios = Some(all_ratios);
        // Clear stale caches — the isotope set changed.
        self.precomputed_cross_sections = None;
        self.precomputed_work_cross_sections = None;
        self.precomputed_base_xs = None;
        // Clear the cubature plan too: atoms are σ-coordinates in
        // ℝ^k and `k` / σ-stack change when groups are reconfigured,
        // so a stale plan would silently produce wrong forward /
        // Jacobian values for the new isotope set even if
        // `cubature_eligible` accepts the same k + grid (Codex round 1
        // P2 on PR #480).
        self.precomputed_sparse_cubature_plan = None;
        self.precomputed_sparse_scalar_plan = None;
        Ok(self)
    }

    // ── Accessors ──

    /// Caller-attached sparse empirical cubature plan, if any.
    /// `spatial_map_typed` reads this so a pre-existing plan is
    /// preserved instead of being clobbered by the local rebuild
    /// pathway (Codex round-5 P3 on PR #480).
    pub fn precomputed_sparse_cubature_plan(
        &self,
    ) -> Option<&Arc<nereids_physics::surrogate::SparseEmpiricalCubaturePlan>> {
        self.precomputed_sparse_cubature_plan.as_ref()
    }

    /// Caller-attached scalar (k = 1) surrogate plan, if any.  Epic
    /// #472, PR #475.
    pub fn precomputed_sparse_scalar_plan(
        &self,
    ) -> Option<&Arc<nereids_physics::surrogate::ScalarSurrogatePlan>> {
        self.precomputed_sparse_scalar_plan.as_ref()
    }

    pub fn energies(&self) -> &[f64] {
        &self.energies
    }
    pub fn resonance_data(&self) -> &[ResonanceData] {
        &self.resonance_data
    }
    pub fn isotope_names(&self) -> &[String] {
        &self.isotope_names
    }
    pub fn temperature_k(&self) -> f64 {
        self.temperature_k
    }
    pub fn resolution(&self) -> Option<&ResolutionFunction> {
        self.resolution.as_ref()
    }
    pub fn initial_densities(&self) -> &[f64] {
        &self.initial_densities
    }
    pub fn solver(&self) -> &SolverConfig {
        &self.solver
    }
    pub fn fit_temperature(&self) -> bool {
        self.fit_temperature
    }
    pub fn transmission_background(&self) -> Option<&BackgroundConfig> {
        self.transmission_background.as_ref()
    }
    pub fn counts_background(&self) -> Option<&CountsBackgroundConfig> {
        self.counts_background.as_ref()
    }
    /// Counts-KL polish override (see [`Self::with_counts_enable_polish`]).
    pub fn counts_enable_polish(&self) -> Option<bool> {
        self.counts_enable_polish
    }
    /// Whether SAMMY TZERO energy-scale calibration is enabled
    /// (see [`Self::with_energy_scale`]).
    pub fn fit_energy_scale(&self) -> bool {
        self.fit_energy_scale
    }
    /// User-specified fit-energy-range restriction (SAMMY EMIN/EMAX
    /// equivalent), or `None` for full-grid fitting.
    /// See [`Self::with_fit_energy_range`].
    pub fn fit_energy_range(&self) -> Option<(f64, f64)> {
        self.fit_energy_range
    }
    pub fn precomputed_cross_sections(&self) -> Option<&Arc<Vec<Vec<f64>>>> {
        self.precomputed_cross_sections.as_ref()
    }
    /// Number of density parameters (one per group or per isotope).
    pub fn n_density_params(&self) -> usize {
        self.n_density_params.unwrap_or(self.resonance_data.len())
    }

    /// Resolve `SolverConfig::Auto` into a concrete solver for the given input.
    pub(crate) fn effective_solver(&self, input: &InputData) -> SolverConfig {
        match &self.solver {
            SolverConfig::Auto => {
                if input.is_counts() {
                    SolverConfig::PoissonKL(PoissonConfig::default())
                } else {
                    SolverConfig::LevenbergMarquardt(LmConfig::default())
                }
            }
            other => other.clone(),
        }
    }
}

/// Fit a single spectrum using the typed input data API.
///
/// Dispatches to the correct fitting engine based on the `InputData` variant
/// and solver configuration:
///
/// | Input | Solver | Path |
/// |-------|--------|------|
/// | Transmission | LM | LM chi-squared with optional SAMMY background |
/// | Transmission | KL | Poisson NLL on transmission with optional background |
/// | Counts | KL | Poisson NLL on raw counts (statistically optimal) |
/// | Counts | LM | Convert to T internally and route to LM |
/// | CountsWithNuisance | KL | Direct Poisson with user-supplied nuisance |
pub fn fit_spectrum_typed(
    input: &InputData,
    config: &UnifiedFitConfig,
) -> Result<SpectrumFitResult, PipelineError> {
    let n_e = config.energies().len();

    // Validate temperature when fitting is requested
    if config.fit_temperature && config.temperature_k < 1.0 {
        return Err(PipelineError::InvalidParameter(format!(
            "temperature must be >= 1.0 K when fit_temperature is true, got {}",
            config.temperature_k,
        )));
    }

    // Validate input length matches energy grid
    if input.n_energies() != n_e {
        return Err(PipelineError::ShapeMismatch(format!(
            "input data has {} energy bins but config.energies has {}",
            input.n_energies(),
            n_e,
        )));
    }

    // Validate auxiliary array lengths match the primary data
    match input {
        InputData::Transmission {
            transmission,
            uncertainty,
        } => {
            if uncertainty.len() != transmission.len() {
                return Err(PipelineError::ShapeMismatch(format!(
                    "uncertainty length {} != transmission length {}",
                    uncertainty.len(),
                    transmission.len(),
                )));
            }
        }
        InputData::Counts {
            sample_counts,
            open_beam_counts,
        } => {
            if open_beam_counts.len() != sample_counts.len() {
                return Err(PipelineError::ShapeMismatch(format!(
                    "open_beam_counts length {} != sample_counts length {}",
                    open_beam_counts.len(),
                    sample_counts.len(),
                )));
            }
        }
        InputData::CountsWithNuisance {
            sample_counts,
            flux,
            background,
        } => {
            if flux.len() != sample_counts.len() {
                return Err(PipelineError::ShapeMismatch(format!(
                    "flux length {} != sample_counts length {}",
                    flux.len(),
                    sample_counts.len(),
                )));
            }
            if background.len() != sample_counts.len() {
                return Err(PipelineError::ShapeMismatch(format!(
                    "background length {} != sample_counts length {}",
                    background.len(),
                    sample_counts.len(),
                )));
            }
        }
    }

    // Reject a malformed caller-supplied precomputed cross-section stack here,
    // before it reaches the forward-model builders (which would otherwise panic
    // on `xs[0].len()` / an over-long row, or silently mis-fit).
    validate_precomputed_cross_sections(config)?;

    let effective_solver = config.effective_solver(input);

    match (input, &effective_solver) {
        // ── Transmission + LM: the well-tested path ──
        (
            InputData::Transmission {
                transmission,
                uncertainty,
            },
            SolverConfig::LevenbergMarquardt(lm_cfg),
        ) => fit_transmission_lm(transmission, uncertainty, config, lm_cfg),

        // ── Transmission + KL: Poisson NLL on transmission values ──
        (
            InputData::Transmission {
                transmission,
                uncertainty,
            },
            SolverConfig::PoissonKL(poisson_cfg),
        ) => fit_transmission_poisson(transmission, uncertainty, config, poisson_cfg),

        // ── Counts + KL: joint-Poisson profile-binomial-deviance path ──
        //
        // The counts-KL solver is now the joint-Poisson fitter validated in
        // memo 35 §P1/§P2 and memo 38.  Uses the explicit `c = Q_s/Q_ob` from
        // `CountsBackgroundConfig::c` and reports `D/(n − k)` as the primary
        // GOF.  Detector-space counts background `B_det` is assumed zero
        // here; the `CountsWithNuisance` arm lets callers supply a
        // detector-bg spectrum.
        (
            InputData::Counts {
                sample_counts,
                open_beam_counts,
            },
            SolverConfig::PoissonKL(poisson_cfg),
        ) => {
            let bg = vec![0.0f64; n_e];
            fit_counts_joint_poisson(
                sample_counts,
                open_beam_counts,
                &bg,
                config,
                &poisson_to_joint_poisson_config(poisson_cfg, config),
            )
        }

        // ── CountsWithNuisance + KL: user-supplied nuisance ──
        (
            InputData::CountsWithNuisance {
                sample_counts,
                flux,
                background,
            },
            SolverConfig::PoissonKL(poisson_cfg),
        ) => fit_counts_joint_poisson(
            sample_counts,
            flux,
            background,
            config,
            &poisson_to_joint_poisson_config(poisson_cfg, config),
        ),

        // ── Counts + LM: convert to transmission (approximate path) ──
        //
        // This is NOT a native counts-domain LM engine.  Counts are divided
        // (sample/OB) to produce transmission, with σ ≈ √max(sample,1)/OB
        // as a simplified Poisson-to-Gaussian conversion.  Poisson structure
        // is lost.  For statistically correct low-count fitting, use the
        // Poisson KL solver (`solver="kl"` or `SolverConfig::Auto`), which
        // now routes to the joint-Poisson path per memo 35 §P1.
        (
            InputData::Counts {
                sample_counts,
                open_beam_counts,
            },
            SolverConfig::LevenbergMarquardt(lm_cfg),
        ) => {
            let (transmission, uncertainty) =
                counts_to_transmission(sample_counts, open_beam_counts);
            fit_transmission_lm(&transmission, &uncertainty, config, lm_cfg)
        }

        // ── CountsWithNuisance + LM: not meaningful ──
        (InputData::CountsWithNuisance { .. }, SolverConfig::LevenbergMarquardt(_)) => {
            Err(PipelineError::InvalidParameter(
                "CountsWithNuisance requires a counts-domain solver (LM cannot use nuisance parameters)"
                    .into(),
            ))
        }

        // Auto should be resolved by effective_solver
        (_, SolverConfig::Auto) => unreachable!("Auto should be resolved before dispatch"),
    }
}

/// Translate the user-facing `PoissonConfig` (payload of `SolverConfig::PoissonKL`)
/// into the internal `JointPoissonFitConfig` required by `joint_poisson_fit`.
///
/// Copies `max_iter` (the only field both structures meaningfully share) and
/// applies any spatial-level polish override carried on the `UnifiedFitConfig`.
fn poisson_to_joint_poisson_config(
    poisson_cfg: &PoissonConfig,
    config: &UnifiedFitConfig,
) -> JointPoissonFitConfig {
    let mut jp_cfg = JointPoissonFitConfig {
        max_iter: poisson_cfg.max_iter,
        ..Default::default()
    };
    if let Some(override_val) = config.counts_enable_polish() {
        jp_cfg.enable_polish = override_val;
    }
    jp_cfg
}

/// Convert counts to transmission: T = sample/open_beam, σ = √(max(sample,1))/open_beam.
///
/// Zero-count bins (sample == 0) get σ = 1e10 so the fitter effectively ignores them.
/// Near-zero open beam bins use a floor of 1e-10 to avoid division by zero.
/// Convert raw counts to transmission with approximate Poisson uncertainty.
///
/// This is a simplified conversion for the Counts+LM fallback path.
/// The uncertainty σ ≈ √max(sample,1)/OB is a Gaussian approximation of
/// Poisson statistics, valid when counts are high (≥ ~20).  At low counts,
/// this overestimates confidence relative to the Poisson KL solver.
///
/// Zero-count and zero-OB bins are marked with sentinel uncertainties
/// (1e10 and 1e30 respectively) so the LM solver effectively ignores them.
fn counts_to_transmission(sample: &[f64], open_beam: &[f64]) -> (Vec<f64>, Vec<f64>) {
    let transmission: Vec<f64> = sample
        .iter()
        .zip(open_beam.iter())
        .map(|(&s, &ob)| if ob > 0.0 { s / ob } else { 0.0 })
        .collect();
    let uncertainty: Vec<f64> = sample
        .iter()
        .zip(open_beam.iter())
        .map(|(&s, &ob)| {
            if ob <= 0.0 {
                // No open beam signal — treat as dead bin
                1e30
            } else if s <= 0.0 {
                // Zero sample counts — large σ so the fitter ignores this bin
                1e10
            } else {
                s.max(1.0).sqrt() / ob.max(1e-10)
            }
        })
        .collect();
    (transmission, uncertainty)
}

/// Transmission + LM path.
fn fit_transmission_lm(
    measured_t: &[f64],
    sigma: &[f64],
    config: &UnifiedFitConfig,
    lm_config: &LmConfig,
) -> Result<SpectrumFitResult, PipelineError> {
    let n_density_params = config.n_density_params();

    // Build parameter vector
    let mut param_vec = build_density_params(config);

    // Guard: energy-scale + temperature fitting is not yet supported.
    // The EnergyScaleTransmissionModel does not wire the temperature parameter.
    if config.fit_energy_scale && config.fit_temperature {
        return Err(PipelineError::InvalidParameter(
            "fit_energy_scale and fit_temperature cannot both be true: \
             EnergyScaleTransmissionModel does not support temperature fitting yet"
                .into(),
        ));
    }

    let _temperature_index = append_temperature_param(&mut param_vec, config);
    let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);

    // Issue #608: seed (t0, L_scale) via resonance peak-matching so the
    // production cold start lands in the global-min basin of the sharply
    // non-convex calibration χ² (the true-σ model's basins are razor-thin).
    seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, measured_t, config);

    // Background parameters (rejects partial BackD/BackF; see
    // validate_transmission_background docstring).
    if let Some(bg) = config.transmission_background.as_ref() {
        validate_transmission_background(bg)?;
    }
    let bg_indices = config
        .transmission_background
        .as_ref()
        .map(|bg| append_background_params(&mut param_vec, bg));

    let mut params = ParameterSet::new(param_vec);
    let mut lm_cfg = lm_config.clone();
    lm_cfg.compute_covariance = config.compute_covariance;

    // Build model — use EnergyScaleTransmissionModel when energy-scale is enabled
    let model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
        build_energy_scale_transmission_model(config, t0_idx, ls_idx)?
    } else {
        build_transmission_model(config, n_density_params, _temperature_index)?
    };

    // Build the per-bin active mask (SAMMY EMIN/EMAX-equivalent fit-energy
    // -range restriction).  `None` when no range is configured — the
    // LM core treats that as "all bins active".
    let active_mask = nereids_fitting::active_mask::build_active_mask(
        config.energies(),
        config.fit_energy_range(),
    );

    // When a fit-energy-range is configured, reject the call early if
    // the user's `[E_min, E_max]` selects fewer active bins than the
    // dispatch can solve.  The LM core's `n_active < n_free` check
    // would also catch the underdetermined case on the main path, but
    // a dispatcher-level rejection gives the user a clear "range too
    // narrow" error instead of a confusing non-converged result, and
    // adds the `n_active < 2` numerical-stability floor so that even
    // an `n_free == 1` fit gets at least one degree of freedom.  See
    // [`required_active_bins`] for the combined `max(2, n_free)`
    // requirement.
    if let Some((e_min, e_max)) = config.fit_energy_range() {
        let n_active = nereids_fitting::active_mask::active_count(
            active_mask.as_deref(),
            config.energies().len(),
        );
        let required = required_active_bins(config);
        if n_active < required {
            return Err(PipelineError::InvalidParameter(format!(
                "fit_energy_range [{e_min}, {e_max}] eV selects {n_active} active bin(s) \
                 on the configured energy grid; at least {required} active bin(s) are \
                 required for LM transmission fitting with {n_free} free parameter(s) \
                 (underdetermined when n_active < n_free)",
                n_free = count_free_params(config),
            )));
        }
    }

    // Dispatch with optional background wrapping
    let result = if let Some(bi) = bg_indices {
        let wrapped = if let (Some(di), Some(fi)) = (bi.back_d, bi.back_f) {
            NormalizedTransmissionModel::new_with_exponential(
                &*model,
                config.energies(),
                bi.anorm,
                bi.back_a,
                bi.back_b,
                bi.back_c,
                di,
                fi,
            )
        } else {
            NormalizedTransmissionModel::new(
                &*model,
                config.energies(),
                bi.anorm,
                bi.back_a,
                bi.back_b,
                bi.back_c,
            )
        };
        lm::levenberg_marquardt_with_mask(
            &wrapped,
            measured_t,
            sigma,
            &mut params,
            &lm_cfg,
            active_mask.as_deref(),
        )?
    } else {
        lm::levenberg_marquardt_with_mask(
            &*model,
            measured_t,
            sigma,
            &mut params,
            &lm_cfg,
            active_mask.as_deref(),
        )?
    };

    let mut sr = extract_result(config, &result, n_density_params, bg_indices)?;

    // Populate energy-scale results if fitted.
    if let Some((t0_idx, ls_idx)) = energy_scale_indices {
        sr.t0_us = Some(result.params[t0_idx]);
        sr.l_scale = Some(result.params[ls_idx]);
    }

    Ok(sr)
}

/// Transmission + Poisson KL path.
///
/// Uses the same model architecture as the LM path:
/// - `EnergyScaleTransmissionModel` when energy-scale fitting is enabled
/// - `NormalizedTransmissionModel` for SAMMY-style background (Anorm + BackA/B/C)
/// - Poisson NLL handles negative model predictions via smooth extrapolation
fn fit_transmission_poisson(
    measured_t: &[f64],
    sigma: &[f64],
    config: &UnifiedFitConfig,
    poisson_cfg: &PoissonConfig,
) -> Result<SpectrumFitResult, PipelineError> {
    // SAMMY EMIN/EMAX-equivalent fit-energy-range (#514): the legacy
    // transmission-domain `poisson_fit` does not honour an active mask,
    // so silently routing the data here when the user set a fit range
    // would bias the fit by including out-of-range bins (margin region)
    // in the cost function.  Hard-reject up-front rather than producing
    // a misleading "successful" fit.  Joint-Poisson (counts) and LM
    // transmission both honour the mask correctly.
    if config.fit_energy_range().is_some() {
        return Err(PipelineError::InvalidParameter(
            "fit_energy_range is not supported for the transmission + \
             Poisson-KL solver path. Use joint-Poisson (provide sample + \
             open-beam counts) or switch to the LM transmission solver."
                .into(),
        ));
    }

    let mut poisson_cfg = poisson_cfg.clone();
    poisson_cfg.compute_covariance = config.compute_covariance;
    let poisson_cfg = &poisson_cfg;

    let n_density_params = config.n_density_params();
    let mut param_vec = build_density_params(config);

    if config.fit_temperature && config.fit_energy_scale {
        return Err(PipelineError::InvalidParameter(
            "fit_energy_scale and fit_temperature cannot both be true: \
             EnergyScaleTransmissionModel does not support temperature fitting yet"
                .into(),
        ));
    }
    let temperature_index = append_temperature_param(&mut param_vec, config);
    let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);

    // Issue #608: peak-match seed for (t0, L_scale) — see fit_transmission_lm.
    seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, measured_t, config);

    // Background parameters — use same SAMMY-style model as LM, with the
    // same partial-BackD/BackF rejection.
    if let Some(bg) = config.transmission_background.as_ref() {
        validate_transmission_background(bg)?;
    }
    let bg_indices = config
        .transmission_background
        .as_ref()
        .map(|bg| append_background_params(&mut param_vec, bg));

    let mut params = ParameterSet::new(param_vec);

    // Build inner model (energy-scale or precomputed)
    let model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
        build_energy_scale_transmission_model(config, t0_idx, ls_idx)?
    } else {
        build_transmission_model(config, n_density_params, temperature_index)?
    };

    // Wrap with NormalizedTransmissionModel for background (same as LM)
    let result = if let Some(bi) = bg_indices {
        let wrapped = if let (Some(di), Some(fi)) = (bi.back_d, bi.back_f) {
            NormalizedTransmissionModel::new_with_exponential(
                &*model,
                config.energies(),
                bi.anorm,
                bi.back_a,
                bi.back_b,
                bi.back_c,
                di,
                fi,
            )
        } else {
            NormalizedTransmissionModel::new(
                &*model,
                config.energies(),
                bi.anorm,
                bi.back_a,
                bi.back_b,
                bi.back_c,
            )
        };
        let pr = poisson::poisson_fit(&wrapped, measured_t, &mut params, poisson_cfg)?;
        poisson_to_lm_result(&wrapped, measured_t, sigma, &pr, &params)
    } else {
        let pr = poisson::poisson_fit(&*model, measured_t, &mut params, poisson_cfg)?;
        poisson_to_lm_result(&*model, measured_t, sigma, &pr, &params)
    }?;

    let mut sr = extract_result(config, &result, n_density_params, bg_indices)?;

    // Populate temperature
    if let Some(idx) = temperature_index {
        sr.temperature_k = Some(result.params[idx]);
        sr.temperature_k_unc = temperature_index.and_then(|i| {
            result
                .uncertainties
                .as_ref()
                .and_then(|u| u.get(i).copied())
        });
    }

    // Populate energy-scale results
    if let Some((t0_idx, ls_idx)) = energy_scale_indices {
        sr.t0_us = Some(result.params[t0_idx]);
        sr.l_scale = Some(result.params[ls_idx]);
    }

    Ok(sr)
}

/// Joint-Poisson counts-path fitter (memo 35 §P1/§P2).
///
/// Builds a pure transmission `FitModel` (density + optional temperature +
/// optional energy-scale) and feeds it to [`joint_poisson::joint_poisson_fit`]
/// together with explicit `(O, S, c)`.  Returns a [`SpectrumFitResult`] with
/// `deviance_per_dof = Some(...)` as the primary GOF (memo 35 §P1.2).
/// `reduced_chi_squared` is set to the same value so GUI consumers that
/// still read the legacy field see a deviance-based metric.
///
/// Current scope (P1 + P2, including P2.2): `fit_alpha_1`, `fit_alpha_2`,
/// and non-zero `detector_background` remain rejected (`λ̂` absorbs the
/// global flux scale, `B_det` / alpha_2 wiring is memo 35 §P3.2 deferred).
/// `transmission_background` with `A_n` + `B_A` / `B_B` / `B_C` is
/// supported as of P2.2, subject to the operational rule that `B_A` must
/// be enabled if any of `B_A` / `B_B` / `B_C` is enabled (memo 35 §P2.2,
/// EG2 S2 C_An shows A_n alone cannot absorb a constant offset — density
/// bias −23%).  Exponential-tail terms `BackD` / `BackF` are rejected
/// (memo 35 §P4-deferred).
fn fit_counts_joint_poisson(
    sample_counts: &[f64],
    flux: &[f64],
    detector_background: &[f64],
    config: &UnifiedFitConfig,
    jp_cfg: &JointPoissonFitConfig,
) -> Result<SpectrumFitResult, PipelineError> {
    // ── Compatibility gates (memo 35 §P3 items are out of scope here) ──
    if let Some(bg) = config.counts_background()
        && (bg.fit_alpha_1 || bg.fit_alpha_2)
    {
        return Err(PipelineError::InvalidParameter(
            "joint-Poisson solver does not support fit_alpha_1/fit_alpha_2: \
             the profile lambda-hat absorbs the global flux scale (alpha_1 redundant); \
             alpha_2 / B_det wiring is deferred to memo 35 §P3."
                .into(),
        ));
    }
    if detector_background.iter().any(|&v| v.abs() > 1e-12) {
        return Err(PipelineError::InvalidParameter(
            "joint-Poisson solver with non-zero detector_background is not yet supported \
             (B_det wiring deferred to memo 35 §P3.2)."
                .into(),
        ));
    }

    // ── §P2.2 operational rule: B_A required if any additive term enabled ──
    if let Some(bg) = config.transmission_background.as_ref() {
        if bg.fit_back_d || bg.fit_back_f {
            return Err(PipelineError::InvalidParameter(
                "joint-Poisson solver does not support the BackD/BackF exponential \
                 tail (memo 35 §P4-deferred)."
                    .into(),
            ));
        }
        if (bg.fit_back_b || bg.fit_back_c) && !bg.fit_back_a {
            return Err(PipelineError::InvalidParameter(
                "joint-Poisson transmission_background: B_A (fit_back_a) must be \
                 enabled whenever any of B_B / B_C is enabled (memo 35 §P2.2 — \
                 A_n alone cannot absorb a constant offset; EG2 S2 C_An → −23% \
                 density bias)."
                    .into(),
            ));
        }
    }

    let c = config.counts_background().map(|b| b.c).unwrap_or(1.0);
    if !(c.is_finite() && c > 0.0) {
        return Err(PipelineError::InvalidParameter(format!(
            "joint-Poisson solver requires finite c > 0 in CountsBackgroundConfig, got {c}",
        )));
    }

    // ── Build parameter vector (density → temperature → energy-scale →
    // transmission background), using the shared
    // append_temperature_param / append_energy_scale_params /
    // append_background_params helpers.
    let n_density_params = config.n_density_params();
    let mut param_vec = build_density_params(config);

    if config.fit_temperature && config.fit_energy_scale {
        return Err(PipelineError::InvalidParameter(
            "fit_energy_scale and fit_temperature cannot both be true: \
             EnergyScaleTransmissionModel does not support temperature fitting yet"
                .into(),
        ));
    }
    let temperature_index = append_temperature_param(&mut param_vec, config);
    let energy_scale_indices = append_energy_scale_params(&mut param_vec, config);

    // Issue #608: peak-match seed for (t0, L_scale).  The KL/counts path fits
    // the same sharply non-convex calibration landscape as the LM path, so it
    // needs the same seed; run it on a (sample − bg)/(flux − bg) transmission
    // proxy (the resonance-dip positions are all the seed needs).
    if energy_scale_indices.is_some() {
        let t_proxy: Vec<f64> = sample_counts
            .iter()
            .zip(flux.iter())
            .zip(detector_background.iter())
            .map(|((&s, &f), &b)| {
                let den = f - b;
                if den > 0.0 {
                    ((s - b).max(0.0) / den).min(2.0)
                } else {
                    1.0
                }
            })
            .collect();
        seed_energy_scale_in_params(&mut param_vec, energy_scale_indices, &t_proxy, config);
    }

    // ── Transmission background (A_n + B_A/B/C) parameters, P2.2 ──
    // Use the same SAMMY-style param block as the LM transmission path.
    // If BackD/BackF were enabled, we would have already errored out above.
    let bg_indices = config
        .transmission_background
        .as_ref()
        .map(|bg| append_background_params(&mut param_vec, bg));

    let mut params = ParameterSet::new(param_vec);

    // ── Build pure transmission model ──
    let t_model: Box<dyn FitModel> = if let Some((t0_idx, ls_idx)) = energy_scale_indices {
        build_energy_scale_transmission_model(config, t0_idx, ls_idx)?
    } else {
        build_transmission_model(config, n_density_params, temperature_index)?
    };

    // ── Wrap with NormalizedTransmissionModel if bg is active (P2.2) ──
    // The wrapper adds `T_out = A_n · T_inner + B_A + B_B/√E + B_C·√E`,
    // exactly matching the SAMMY form used by the LM transmission path.
    // Its analytical Jacobian chains through the inner model correctly,
    // so JointPoissonObjective picks up gradients for both density and
    // background parameters without further wiring.

    // Build the per-bin active mask (SAMMY EMIN/EMAX-equivalent fit-energy
    // -range restriction).  `None` when no range is configured — the
    // JP objective treats that as "all bins active".
    let active_mask = nereids_fitting::active_mask::build_active_mask(
        config.energies(),
        config.fit_energy_range(),
    );
    let active_mask_slice = active_mask.as_deref();

    // Reject early when the configured range selects fewer active
    // bins than the joint-Poisson dispatch can solve.  `joint_poisson_fit`
    // already rejects the `n_active < n_free` case (and the
    // `n_active == 0` early-return), but a dispatcher-level rejection
    // gives users a clear "range too narrow" error instead of a
    // non-converged JP result, and adds the `n_active < 2` numerical-
    // stability floor so even an `n_free == 1` fit gets at least one
    // degree of freedom.  See [`required_active_bins`] for the
    // combined `max(2, n_free)` requirement.
    if let Some((e_min, e_max)) = config.fit_energy_range() {
        let n_active =
            nereids_fitting::active_mask::active_count(active_mask_slice, config.energies().len());
        let required = required_active_bins(config);
        if n_active < required {
            return Err(PipelineError::InvalidParameter(format!(
                "fit_energy_range [{e_min}, {e_max}] eV selects {n_active} active bin(s) \
                 on the configured energy grid; at least {required} active bin(s) are \
                 required for joint-Poisson fitting with {n_free} free parameter(s) \
                 (underdetermined when n_active < n_free)",
                n_free = count_free_params(config),
            )));
        }
    }

    let result;
    if let Some(bi) = bg_indices {
        let wrapped = NormalizedTransmissionModel::new(
            &*t_model,
            config.energies(),
            bi.anorm,
            bi.back_a,
            bi.back_b,
            bi.back_c,
        );
        let objective = JointPoissonObjective {
            model: &wrapped,
            o: flux,
            s: sample_counts,
            c,
            active_mask: active_mask_slice,
        };
        let mut cfg = jp_cfg.clone();
        cfg.compute_covariance = config.compute_covariance;
        result = joint_poisson::joint_poisson_fit(&objective, &mut params, &cfg).map_err(|e| {
            PipelineError::InvalidParameter(format!("joint-Poisson fit failed: {e}"))
        })?;
    } else {
        let objective = JointPoissonObjective {
            model: &*t_model,
            o: flux,
            s: sample_counts,
            c,
            active_mask: active_mask_slice,
        };
        let mut cfg = jp_cfg.clone();
        cfg.compute_covariance = config.compute_covariance;
        result = joint_poisson::joint_poisson_fit(&objective, &mut params, &cfg).map_err(|e| {
            PipelineError::InvalidParameter(format!("joint-Poisson fit failed: {e}"))
        })?;
    }

    // ── Extract fitted quantities ──
    let densities: Vec<f64> = (0..n_density_params).map(|i| result.params[i]).collect();

    let (uncertainties, temperature_k_unc) = if let Some(ref unc_all) = result.uncertainties {
        let dens_unc: Vec<f64> = (0..n_density_params)
            .map(|i| *unc_all.get(i).unwrap_or(&f64::NAN))
            .collect();
        let t_unc = temperature_index.and_then(|idx| {
            params
                .free_indices()
                .iter()
                .position(|&fi| fi == idx)
                .and_then(|pos| unc_all.get(pos).copied())
        });
        (Some(dens_unc), t_unc)
    } else {
        (None, None)
    };
    let fitted_temp = temperature_index.map(|idx| result.params[idx]);

    // Convergence signal per memo 35 §P2.3: the deviance value is the
    // acceptance criterion, but we expose a boolean to preserve the
    // existing SpectrumFitResult shape.  Report True when EITHER stage
    // self-flagged convergence (whichever accepts).
    let converged = result.gn_converged || result.polish_converged;

    // ── Background parameter readout ──
    // When bg is active, read A_n / B_A / B_B / B_C from the fitted
    // parameter vector at their registered indices.  When bg is absent,
    // use the memo 35 §P1 convention: A_n = 1 (subsumed into λ̂), bg = 0.
    let (anorm_out, bg_abc_out) = if let Some(bi) = bg_indices {
        (
            result.params[bi.anorm],
            [
                result.params[bi.back_a],
                result.params[bi.back_b],
                result.params[bi.back_c],
            ],
        )
    } else {
        (1.0, [0.0, 0.0, 0.0])
    };

    Ok(SpectrumFitResult {
        densities,
        uncertainties,
        // Back-compat bridge: reduced_chi_squared carries D/(n−k) for the
        // joint-Poisson path.  Memo 35 §P1.2 — Pearson χ² is secondary.
        reduced_chi_squared: result.deviance_per_dof,
        converged,
        iterations: result.gn_iterations + result.polish_iterations,
        temperature_k: fitted_temp,
        temperature_k_unc,
        anorm: anorm_out,
        background: bg_abc_out,
        // Joint-Poisson never fits the exponential tail.  `None`
        // signals "tail not fit" to downstream consumers (GUI overlay
        // curve drops the exponential term; PyO3 conversion passes
        // through to Python as `None`).
        back_d: None,
        back_f: None,
        t0_us: energy_scale_indices.map(|(t0_idx, _)| result.params[t0_idx]),
        l_scale: energy_scale_indices.map(|(_, ls_idx)| result.params[ls_idx]),
        deviance_per_dof: Some(result.deviance_per_dof),
    })
}

// ── Shared helpers for fit_spectrum_typed ──

fn build_density_params(config: &UnifiedFitConfig) -> Vec<FitParameter> {
    config
        .initial_densities
        .iter()
        .enumerate()
        .map(|(i, &d)| {
            FitParameter::non_negative(
                config
                    .isotope_names
                    .get(i)
                    .cloned()
                    .unwrap_or_else(|| format!("isotope_{i}")),
                d,
            )
        })
        .collect()
}

/// Append a temperature parameter to the fit vector if
/// `config.fit_temperature` is `true`.  Returns the parameter index.
///
/// Bounds [1.0, 5000.0] K match the transmission / counts fit paths.
/// Temperature unit is Kelvin; the initial value is `config.temperature_k`.
fn append_temperature_param(
    param_vec: &mut Vec<FitParameter>,
    config: &UnifiedFitConfig,
) -> Option<usize> {
    if !config.fit_temperature {
        return None;
    }
    let idx = param_vec.len();
    param_vec.push(FitParameter {
        name: "temperature_k".into(),
        value: config.temperature_k,
        lower: 1.0,
        upper: 5000.0,
        fixed: false,
    });
    Some(idx)
}

/// Detect significant transmission dips (resonance signatures) in a measured
/// spectrum: returns `(energy_eV, depth)` for each local minimum whose depth
/// below the no-absorption baseline is a meaningful fraction of the deepest
/// dip.  Light 3-point smoothing suppresses single-bin noise.  Used to seed the
/// energy-scale calibration (issue #608).
fn detect_transmission_dips(measured: &[f64], energies: &[f64]) -> Vec<(f64, f64)> {
    let n = measured.len();
    if n < 5 || energies.len() != n {
        return Vec::new();
    }
    let mut sm = measured.to_vec();
    for i in 1..n - 1 {
        sm[i] = (measured[i - 1] + measured[i] + measured[i + 1]) / 3.0;
    }
    // No-absorption baseline ≈ a high percentile of the (smoothed) signal.
    let mut sorted = sm.clone();
    sorted.sort_by(f64::total_cmp);
    let baseline = sorted[((0.9 * (n as f64 - 1.0)).round() as usize).min(n - 1)];
    let mut dips: Vec<(f64, f64)> = Vec::new();
    let mut max_depth = 0.0f64;
    for i in 1..n - 1 {
        if sm[i] < sm[i - 1] && sm[i] <= sm[i + 1] {
            let depth = baseline - sm[i];
            if depth > 0.0 {
                dips.push((energies[i], depth));
                max_depth = max_depth.max(depth);
            }
        }
    }
    // Keep dips that are a meaningful fraction of the deepest — filters noise
    // wiggles while retaining genuine (even weak) resonance dips.
    dips.retain(|&(_, d)| d >= 0.2 * max_depth);
    dips
}

/// Seed `(t0, L_scale)` for the energy-scale fit by resonance peak-matching.
///
/// A TOF calibration is linear: `tof_measured = t0 + L_scale · tof_nominal`.
/// We detect the transmission dips in the measured spectrum, match each to its
/// nearest known resonance energy (under the nominal calibration), and
/// least-squares-fit the matched `(tof_nominal, tof_measured)` pairs for
/// `(t0, L_scale)`.  This is landscape-independent — it seeds the LM into the
/// global-minimum basin of the sharply non-convex post-#608 calibration χ²,
/// which a cold start (t0=0, L_scale=1) or a grid scan cannot reliably find
/// (the physically-exact true-σ model's basins are razor-thin).  The physics is
/// unchanged; this only chooses the LM's starting point.
///
/// Returns `None` (→ caller keeps the configured cold start) when fewer than
/// two distinct resonances can be matched, so it is a safe enhancement: it
/// improves the seed when it can and never degrades the existing behaviour.
fn peak_match_energy_scale_seed(
    measured: &[f64],
    energies: &[f64],
    config: &UnifiedFitConfig,
    flight_path_m: f64,
    t0_bounds: (f64, f64),
    l_scale_bounds: (f64, f64),
) -> Option<(f64, f64)> {
    let n = energies.len();
    if n < 5 || measured.len() != n {
        return None;
    }
    // Resonance centers within the measured energy range (true frame).
    let refs: Vec<&ResonanceData> = config.resonance_data().iter().collect();
    let (e_lo, e_hi) = (
        energies[0].min(energies[n - 1]),
        energies[0].max(energies[n - 1]),
    );
    let res_e: Vec<f64> = nereids_physics::transmission::resonance_center_energies(&refs)
        .into_iter()
        .filter(|&e| e > e_lo && e < e_hi)
        .collect();
    if res_e.len() < 2 {
        return None;
    }
    let dips = detect_transmission_dips(measured, energies);
    if dips.len() < 2 {
        return None;
    }
    // Match each dip to its nearest resonance (nominal calibration ⇒
    // corrected ≈ measured), keeping the deepest dip per resonance.  Reject a
    // dip that does not sit UNAMBIGUOUSLY near a resonance — within half the
    // minimum inter-resonance spacing — so a spurious/mis-matched dip drops out
    // rather than poisoning the linear fit (issue #608 R2).  `res_e` is sorted
    // and de-duplicated.
    //
    // Floor `match_tol` at the energy-grid resolution (Copilot PR #609): a dip
    // is localized to ~one grid step, so a single closely-spaced resonance pair
    // must not collapse the GLOBAL tolerance toward 0 and reject dips at
    // well-separated resonances.  Dedup already stops exact duplicates from
    // zeroing it; the floor additionally guards distinct near-degenerate
    // energies.  Well-separated resonances keep `0.5·min_spacing ≫ grid_res`, so
    // the floor is inert in the common case.
    let min_spacing = res_e
        .windows(2)
        .map(|w| (w[1] - w[0]).abs())
        .fold(f64::INFINITY, f64::min);
    let mut grid_steps: Vec<f64> = energies.windows(2).map(|w| (w[1] - w[0]).abs()).collect();
    grid_steps.sort_by(f64::total_cmp);
    let grid_res = grid_steps.get(grid_steps.len() / 2).copied().unwrap_or(0.0);
    let match_tol = (0.5 * min_spacing).max(grid_res);
    let mut best: Vec<Option<(f64, f64)>> = vec![None; res_e.len()];
    for &(e_dip, depth) in &dips {
        let Some((k, &re)) = res_e
            .iter()
            .enumerate()
            .min_by(|(_, a), (_, b)| (*a - e_dip).abs().total_cmp(&(*b - e_dip).abs()))
        else {
            continue;
        };
        if (re - e_dip).abs() > match_tol {
            continue;
        }
        match best[k] {
            Some((_, d)) if d >= depth => {}
            _ => best[k] = Some((e_dip, depth)),
        }
    }
    // Matched (tof_nominal, tof_measured) pairs.  The common factor
    // `tof_factor · flight_path` must match `EnergyScaleTransmissionModel`'s
    // `corrected_energies` (it cancels in the slope but sets the intercept t0).
    let tof_factor = (0.5 * NEUTRON_MASS_KG / EV_TO_JOULES).sqrt() * 1.0e6;
    let c = tof_factor * flight_path_m;
    let pairs: Vec<(f64, f64)> = res_e
        .iter()
        .zip(best.iter())
        .filter_map(|(&re, b)| b.map(|(e_dip, _)| (c / re.sqrt(), c / e_dip.sqrt())))
        .collect();
    if pairs.len() < 2 {
        return None;
    }
    // Linear least squares: tof_meas = t0 + L_scale · tof_nom.
    let m = pairs.len() as f64;
    let mean_x = pairs.iter().map(|p| p.0).sum::<f64>() / m;
    let mean_y = pairs.iter().map(|p| p.1).sum::<f64>() / m;
    let mut sxx = 0.0f64;
    let mut sxy = 0.0f64;
    for &(x, y) in &pairs {
        sxx += (x - mean_x) * (x - mean_x);
        sxy += (x - mean_x) * (y - mean_y);
    }
    if sxx <= 0.0 {
        return None;
    }
    let l_scale = sxy / sxx;
    let t0 = mean_y - l_scale * mean_x;
    if !t0.is_finite() || !l_scale.is_finite() {
        return None;
    }
    // Reject (→ keep the configured cold start) when the fitted seed falls
    // outside the parameter bounds.  Clamping an out-of-range fit onto a bound
    // would seed the LM worse than its cold start — the opposite of this seed's
    // purpose (issue #608 R2).  An in-range fit is the high-confidence case.
    if t0 < t0_bounds.0
        || t0 > t0_bounds.1
        || l_scale < l_scale_bounds.0
        || l_scale > l_scale_bounds.1
    {
        return None;
    }
    Some((t0, l_scale))
}

/// Apply the peak-matching seed to the `(t0, L_scale)` entries of `param_vec`
/// in place, when energy-scale fitting is enabled (issue #608).  No-op when the
/// seed cannot be computed (→ the configured cold start is kept).
fn seed_energy_scale_in_params(
    param_vec: &mut [FitParameter],
    energy_scale_indices: Option<(usize, usize)>,
    measured_transmission: &[f64],
    config: &UnifiedFitConfig,
) {
    let Some((t0_idx, ls_idx)) = energy_scale_indices else {
        return;
    };
    let t0_b = (param_vec[t0_idx].lower, param_vec[t0_idx].upper);
    let ls_b = (param_vec[ls_idx].lower, param_vec[ls_idx].upper);
    if let Some((t0_seed, ls_seed)) = peak_match_energy_scale_seed(
        measured_transmission,
        config.energies(),
        config,
        config.flight_path_m,
        t0_b,
        ls_b,
    ) {
        param_vec[t0_idx].value = t0_seed;
        param_vec[ls_idx].value = ls_seed;
    }
}

/// Append SAMMY TZERO energy-scale parameters (t_0 and L_scale) when
/// `config.fit_energy_scale` is `true`.  Returns `(t0_idx, l_scale_idx)`.
///
/// Bounds: `t_0 ∈ [-10.0, 10.0] μs` and `L_scale ∈ [0.99, 1.01]`
/// (dimensionless).  Matches the LM / KL transmission paths.
fn append_energy_scale_params(
    param_vec: &mut Vec<FitParameter>,
    config: &UnifiedFitConfig,
) -> Option<(usize, usize)> {
    if !config.fit_energy_scale {
        return None;
    }
    let t0_idx = param_vec.len();
    param_vec.push(FitParameter {
        name: "t0_us".into(),
        value: config.t0_init_us,
        lower: -10.0,
        upper: 10.0,
        fixed: false,
    });
    let ls_idx = param_vec.len();
    param_vec.push(FitParameter {
        name: "l_scale".into(),
        value: config.l_scale_init,
        lower: 0.99,
        upper: 1.01,
        fixed: false,
    });
    Some((t0_idx, ls_idx))
}

/// Validate that the SAMMY-style `BackgroundConfig` is internally
/// consistent for the purposes of the production fit dispatch.
///
/// **BackD / BackF must travel together.**  The
/// `NormalizedTransmissionModel` exponential-tail wrapper takes both
/// indices or neither (`new_with_exponential` requires both, `new` takes
/// neither).  Allowing only one of `fit_back_d` / `fit_back_f` would
/// leave the other parameter registered as "free" but absent from the
/// objective and Jacobian — silently fitting nothing while reporting a
/// misleading initial value back to the caller.  Reject the partial
/// configuration up-front with a clear error.
///
/// Idempotent on valid configs; intended to be called by every
/// production fit entry that takes a `transmission_background`.
pub(crate) fn validate_transmission_background(bg: &BackgroundConfig) -> Result<(), PipelineError> {
    if bg.fit_back_d != bg.fit_back_f {
        return Err(PipelineError::InvalidParameter(format!(
            "transmission_background: fit_back_d ({}) and fit_back_f ({}) \
             must both be true or both be false. The exponential tail \
             wrapper (BackD · exp(−BackF / √E)) requires both parameters \
             together; enabling only one leaves the other registered but \
             unused, silently producing the initial value as the fitted \
             result. Either enable both (to fit the exponential tail) or \
             disable both (4-term wrapper without exponential).",
            bg.fit_back_d, bg.fit_back_f,
        )));
    }
    Ok(())
}

/// Validate a caller-supplied precomputed cross-section stack against the
/// config's grid and isotope/group mapping, BEFORE it reaches the forward-model
/// builders or the per-pixel rayon loop.
///
/// `precomputed_cross_sections` is consumed by `build_transmission_model` /
/// `build_energy_scale_transmission_model`, which index `xs[0].len()` (panics
/// when empty) and `params[density_indices[i]]` for each row, and write
/// `neg_opt[j]` for every `j` in a row (an over-long row writes out of bounds).
/// A shape mismatch there either panics deep in the LM iteration or is swallowed
/// per-pixel as `n_failed`; validating once up front turns it into a typed
/// `ShapeMismatch` (mapped to `PyValueError` at the Python boundary).
///
/// Accepted shapes (matching the builders' collapse logic):
/// * **non-empty** — at least one σ row.
/// * every row has length `energies.len()`.
/// * row count is either `n_density_params` (already group-collapsed / identity)
///   or, when groups are active, `density_indices.len()` (per-member, collapsed
///   downstream).
///
/// When `precomputed_work_cross_sections` is also set (issue #608, Gaussian
/// aux-grid path) its working-grid σ + layout are validated against the same
/// invariants: non-empty, each row length == `work_layout.energies.len()`,
/// finite σ, row count == the data-grid σ row count (same density mapping), and
/// a layout whose `data_indices` length == `energies.len()` with every index in
/// range for the working grid — so `build_transmission_model`'s Beer-Lambert
/// accumulation and `work_layout.extract(..)` cannot write/read out of bounds.
pub(crate) fn validate_precomputed_cross_sections(
    config: &UnifiedFitConfig,
) -> Result<(), PipelineError> {
    let Some(xs) = config.precomputed_cross_sections() else {
        return Ok(());
    };
    if xs.is_empty() {
        return Err(PipelineError::ShapeMismatch(
            "precomputed_cross_sections must not be empty".into(),
        ));
    }

    let n_e = config.energies().len();
    for (i, row) in xs.iter().enumerate() {
        if row.len() != n_e {
            return Err(PipelineError::ShapeMismatch(format!(
                "precomputed_cross_sections row {i} has length {} but config.energies \
                 has {n_e}",
                row.len(),
            )));
        }
        // A correctly-shaped row of NaN / ±∞ σ passes the shape checks but
        // poisons the forward model: every transmission sample picks up the
        // non-finite σ and the LM residual / Fisher matrix becomes NaN, which
        // the per-pixel dispatch then silently swallows as a failed fit.
        // Reject non-finite σ at the boundary (`is_finite()` excludes both NaN
        // and ±∞; a bare order comparison would let NaN through).
        if let Some(j) = row.iter().position(|s| !s.is_finite()) {
            return Err(PipelineError::ShapeMismatch(format!(
                "precomputed_cross_sections row {i} has non-finite σ at energy index {j}: {}",
                row[j],
            )));
        }
    }

    let n_params = config.n_density_params();
    // When groups are active the per-member form (`density_indices.len()` rows)
    // is collapsed to `n_params` rows downstream; accept either.
    let member_rows = config.density_indices.as_ref().map(|di| di.len());
    let row_ok = xs.len() == n_params || member_rows == Some(xs.len());
    if !row_ok {
        let expected = match member_rows {
            Some(m) if m != n_params => format!("{n_params} (collapsed) or {m} (per-member)"),
            _ => format!("{n_params}"),
        };
        return Err(PipelineError::ShapeMismatch(format!(
            "precomputed_cross_sections has {} rows but expected {expected}",
            xs.len(),
        )));
    }

    // Issue #608: the working-grid σ + layout attached via
    // `with_precomputed_work_cross_sections` flows into
    // `build_transmission_model`, where the model applies Beer-Lambert +
    // resolution on `work_layout.energies` and then `work_layout.extract(..)`
    // indexes the broadened spectrum by `work_layout.data_indices`.  An empty
    // σ panics on `xs[0].len()`; a row whose length ≠ the working-grid length
    // writes out of bounds in the Beer-Lambert accumulation; a layout whose
    // `data_indices` length ≠ the data grid, or that indexes past the working
    // grid, panics in `extract`.  Validate the same shape/consistency
    // invariants as the data-grid σ above so a malformed setter call surfaces
    // as a typed `ShapeMismatch` here rather than a per-pixel panic / swallowed
    // failed fit.
    if let Some((work_xs, layout)) = &config.precomputed_work_cross_sections {
        let n_work = layout.energies.len();
        if work_xs.is_empty() {
            return Err(PipelineError::ShapeMismatch(
                "precomputed_work_cross_sections must not be empty".into(),
            ));
        }
        for (i, row) in work_xs.iter().enumerate() {
            if row.len() != n_work {
                return Err(PipelineError::ShapeMismatch(format!(
                    "precomputed_work_cross_sections row {i} has length {} but \
                     work_layout has {n_work} working-grid energies",
                    row.len(),
                )));
            }
            if let Some(j) = row.iter().position(|s| !s.is_finite()) {
                return Err(PipelineError::ShapeMismatch(format!(
                    "precomputed_work_cross_sections row {i} has non-finite σ at \
                     working-grid index {j}: {}",
                    row[j],
                )));
            }
        }
        // Row count must match the data-grid σ row count (same density mapping):
        // both feed the SAME `density_indices` in `build_transmission_model`.
        if work_xs.len() != xs.len() {
            return Err(PipelineError::ShapeMismatch(format!(
                "precomputed_work_cross_sections has {} rows but \
                 precomputed_cross_sections has {} — both index the same density \
                 mapping and must agree",
                work_xs.len(),
                xs.len(),
            )));
        }
        // The layout maps each data energy to a working-grid index.  Its length
        // must equal the data grid, and every index must be in range so
        // `extract` cannot panic.
        if layout.data_indices.len() != n_e {
            return Err(PipelineError::ShapeMismatch(format!(
                "precomputed_work_cross_sections layout maps {} data points but \
                 config.energies has {n_e}",
                layout.data_indices.len(),
            )));
        }
        if let Some(&bad) = layout.data_indices.iter().find(|&&idx| idx >= n_work) {
            return Err(PipelineError::ShapeMismatch(format!(
                "precomputed_work_cross_sections layout index {bad} is out of \
                 range for {n_work} working-grid energies",
            )));
        }
    }
    Ok(())
}

/// Count the number of free parameters the production LM transmission
/// or joint-Poisson (counts-KL) dispatch will register for `config`.
///
/// Mirrors the parameter-vector assembly performed by
/// [`fit_transmission_lm`], [`fit_counts_joint_poisson`] and the
/// shared `build_density_params` / `append_*_param` helpers below:
///
/// * `n_density_params` density slots (always free).
/// * `+1` if [`UnifiedFitConfig::fit_temperature`] is set.
/// * `+2` if [`UnifiedFitConfig::fit_energy_scale`] is set
///   (t_0 and L_scale).
/// * For a transmission_background, the count of `true` flags on
///   `(fit_anorm, fit_back_a, fit_back_b, fit_back_c, fit_back_d,
///   fit_back_f)`.  The joint-Poisson dispatch rejects BackD/BackF
///   up-front, but this helper counts them as free when set so the
///   fail-fast diagnostic ordering (`underdetermined > BackD/BackF
///   reject`) does not flip across paths.
///
/// Used both by [`validate_spatial_fit_preflight`] (whole-map
/// rejection) and by the per-pixel LM / JP dispatchers
/// (single-spectrum rejection) so the `n_active < required` bound
/// stays in lockstep across the spatial / single-spectrum boundary.
/// The research-only `fit_alpha_1` / `fit_alpha_2` flags on
/// `CountsBackgroundConfig` are *not* counted: the joint-Poisson
/// dispatch rejects them up-front and the LM path never wires them.
pub(crate) fn count_free_params(config: &UnifiedFitConfig) -> usize {
    let mut n_free = config.n_density_params();
    if config.fit_temperature {
        n_free += 1;
    }
    if config.fit_energy_scale {
        n_free += 2;
    }
    if let Some(bg) = config.transmission_background.as_ref() {
        n_free += usize::from(bg.fit_anorm);
        n_free += usize::from(bg.fit_back_a);
        n_free += usize::from(bg.fit_back_b);
        n_free += usize::from(bg.fit_back_c);
        n_free += usize::from(bg.fit_back_d);
        n_free += usize::from(bg.fit_back_f);
    }
    n_free
}

/// Minimum number of active bins the LM / joint-Poisson dispatch
/// requires for a fit on `config`.  Encodes two distinct constraints:
///
/// 1. **Underdetermined-system rejection**: a fit with `n_active <
///    n_free` cannot be solved (the Jacobian cannot be full rank).
///    Lifting this check out of the LM / joint-Poisson cores into
///    the dispatcher gives the user an actionable error instead of
///    silent NaN propagation through the spatial map.
/// 2. **Numerical-stability floor**: even with `n_free == 1`, a
///    one-active-bin fit (`n_active == 1`) is exactly determined
///    (`dof == 0`) and gives no curvature for the LM step / no
///    deviance signal for the joint-Poisson reduced GOF.  We keep
///    the previous `n_active < 2` minimum so callers that relied
///    on it (e.g. the spatial-preflight regression tests for
///    narrow `fit_energy_range` windows) continue to receive the
///    same actionable diagnostic.
///
/// The combined requirement is `max(2, count_free_params(config))`.
pub(crate) fn required_active_bins(config: &UnifiedFitConfig) -> usize {
    count_free_params(config).max(2)
}

fn append_background_params(
    param_vec: &mut Vec<FitParameter>,
    bg: &BackgroundConfig,
) -> BackgroundIndices {
    // Anorm bounded to [0.5, 2.0] — physically reasonable normalization range.
    // Previously unbounded [0, ∞), which allowed the fitter to absorb signal
    // into anorm (e.g., anorm=15.9 with density=0.03×true).
    // SAMMY also bounds normalization to a reasonable range.
    let anorm = param_vec.len();
    param_vec.push(if bg.fit_anorm {
        FitParameter {
            name: "anorm".into(),
            value: bg.anorm_init,
            lower: 0.5,
            upper: 2.0,
            fixed: false,
        }
    } else {
        FitParameter::fixed("anorm", bg.anorm_init)
    });
    // Background terms bounded to [-0.5, 0.5].
    // These are small corrections to the transmission baseline.
    // Unbounded background allows the fitter to absorb resonance signal
    // into the background polynomial, producing meaningless densities.
    // SAMMY also constrains background to reasonable ranges.
    let back_a = param_vec.len();
    param_vec.push(if bg.fit_back_a {
        FitParameter {
            name: "back_a".into(),
            value: bg.back_a_init,
            lower: -0.5,
            upper: 0.5,
            fixed: false,
        }
    } else {
        FitParameter::fixed("back_a", bg.back_a_init)
    });
    let back_b = param_vec.len();
    param_vec.push(if bg.fit_back_b {
        FitParameter {
            name: "back_b".into(),
            value: bg.back_b_init,
            lower: -0.5,
            upper: 0.5,
            fixed: false,
        }
    } else {
        FitParameter::fixed("back_b", bg.back_b_init)
    });
    let back_c = param_vec.len();
    param_vec.push(if bg.fit_back_c {
        FitParameter {
            name: "back_c".into(),
            value: bg.back_c_init,
            lower: -0.5,
            upper: 0.5,
            fixed: false,
        }
    } else {
        FitParameter::fixed("back_c", bg.back_c_init)
    });

    // Exponential tail: BackD × exp(−BackF / √E).
    // SAMMY manual Sec III.E.2 — terms 5-6.
    // BackD (amplitude): non-negative, bounded [0, 1].
    // BackF (decay constant in √eV units): non-negative, bounded [0, 100].
    // Note: if BackD_init = 0, the Jacobian column for BackF is identically
    // zero and the optimizer cannot learn BackF.  The default init (0.01)
    // avoids this.
    let back_d = if bg.fit_back_d {
        let idx = param_vec.len();
        param_vec.push(FitParameter {
            name: "back_d".into(),
            value: bg.back_d_init,
            lower: 0.0,
            upper: 1.0,
            fixed: false,
        });
        Some(idx)
    } else {
        None
    };
    let back_f = if bg.fit_back_f {
        let idx = param_vec.len();
        param_vec.push(FitParameter {
            name: "back_f".into(),
            value: bg.back_f_init,
            lower: 0.0,
            upper: 100.0,
            fixed: false,
        });
        Some(idx)
    } else {
        None
    };

    BackgroundIndices {
        anorm,
        back_a,
        back_b,
        back_c,
        back_d,
        back_f,
    }
}

/// Build an [`EnergyScaleTransmissionModel`] for the energy-scale fit branch
/// of `fit_transmission_lm` / `fit_transmission_poisson` /
/// `fit_counts_joint_poisson`.
///
/// Absorbs the 3 byte-identical construction blocks that used to live inline
/// in each fitter (see commit `1b4131f` for the pre-extraction pattern).
/// Stays private — internal pipeline construction detail, not part of the
/// crate's public surface.
///
/// Issue #608: the energy-scale model evaluates the TRUE cross-section at the
/// corrected energies from resonance data (matching `forward_model`) rather
/// than interpolating a precomputed σ grid, so it is constructed from the
/// per-isotope `config.resonance_data`, the density mapping
/// (`config.density_indices` / `config.density_ratios`, defaulting to the
/// identity mapping with ratio 1.0 when no groups are configured), and
/// `config.temperature_k` for Doppler broadening.  Resolution is NOT applied
/// here — the model applies it inside `evaluate()` on the working grid via the
/// wrapped [`InstrumentParams`].  The per-isotope density accumulation
/// (`params[density_indices[i]] * density_ratios[i]`) happens inside the model,
/// so no cross-section pre-collapse is done here.
///
/// `t0_idx` and `ls_idx` are the parameter indices for `t0` and `l_scale`,
/// produced by [`append_energy_scale_params`] at each fitter's setup.
fn build_energy_scale_transmission_model(
    config: &UnifiedFitConfig,
    t0_idx: usize,
    ls_idx: usize,
) -> Result<Box<dyn FitModel>, PipelineError> {
    let instrument = config
        .resolution
        .clone()
        .map(|r| Arc::new(InstrumentParams { resolution: r }));
    // Issue #608: EnergyScale evaluates the TRUE σ at the corrected energies
    // from resonance data (matching forward_model) rather than interpolating a
    // precomputed σ grid, so it takes the per-isotope resonance data + density
    // mapping + temperature.  When no groups are configured, each isotope is its
    // own density parameter (identity mapping, ratio 1.0).
    let n_iso = config.resonance_data.len();
    let density_indices = config
        .density_indices
        .clone()
        .unwrap_or_else(|| (0..n_iso).collect());
    let density_ratios = config
        .density_ratios
        .clone()
        .unwrap_or_else(|| vec![1.0; n_iso]);
    let mut es_model = EnergyScaleTransmissionModel::new(
        Arc::new(config.resonance_data.clone()),
        Arc::new(density_indices),
        Arc::new(density_ratios),
        config.temperature_k,
        config.energies.clone(),
        config.flight_path_m,
        t0_idx,
        ls_idx,
        instrument,
    );
    if let Some(method) = config.tzero_jacobian_method {
        es_model = es_model.with_jacobian_method(method);
    }
    Ok(Box::new(es_model))
}

/// Build the transmission forward model, selecting precomputed or full path.
fn build_transmission_model(
    config: &UnifiedFitConfig,
    n_density_params: usize,
    temperature_index: Option<usize>,
) -> Result<Box<dyn FitModel>, PipelineError> {
    let n_params = config.n_density_params();
    if !config.fit_temperature
        && let Some(xs) = &config.precomputed_cross_sections
    {
        // When groups are active, compute σ_eff per group from member XS.
        // For ungrouped isotopes, this is a no-op (identity mapping, ratio=1.0).
        // Only collapse when XS is per-member (shape matches mapping); if XS is
        // already group-collapsed (len == n_params), this is a clone.
        let collapse_by_groups = |xs: &Arc<Vec<Vec<f64>>>| -> Arc<Vec<Vec<f64>>> {
            if let (Some(di), Some(dr)) = (&config.density_indices, &config.density_ratios)
                && xs.len() == di.len()
                && di.len() == dr.len()
            {
                let n_e = xs[0].len();
                let mut eff = vec![vec![0.0f64; n_e]; n_params];
                for ((&idx, &ratio), member_xs) in di.iter().zip(dr.iter()).zip(xs.iter()) {
                    for (j, &sigma) in member_xs.iter().enumerate() {
                        eff[idx][j] += ratio * sigma;
                    }
                }
                return Arc::new(eff);
            }
            Arc::clone(xs)
        };

        // Issue #608: prefer the WORKING-grid σ + layout when the spatial
        // builder injected it (Gaussian resolution → auxiliary extended grid).
        // The model then applies resolution on the working grid and extracts
        // the data points last.  When absent (tabulated / no resolution) the
        // working grid is the data grid: use the data-grid σ with no layout, so
        // the surrogate fast paths and data-grid `resolution_plan` are
        // unaffected.
        let (effective_xs, work_layout): (
            Arc<Vec<Vec<f64>>>,
            Option<Arc<nereids_physics::transmission::WorkingGridLayout>>,
        ) = match &config.precomputed_work_cross_sections {
            Some((work_xs, layout)) => (collapse_by_groups(work_xs), Some(Arc::clone(layout))),
            None => (collapse_by_groups(xs), None),
        };
        // Issue #442: pass energies + instrument so evaluate() applies
        // resolution after Beer-Lambert on total transmission.
        let instrument = config
            .resolution
            .clone()
            .map(|r| Arc::new(InstrumentParams { resolution: r }));
        let resolution_plan = if instrument.is_some() {
            config.precomputed_resolution_plan.clone()
        } else {
            None
        };
        let sparse_cubature_plan = if instrument.is_some() {
            config.precomputed_sparse_cubature_plan.clone()
        } else {
            None
        };
        let sparse_scalar_plan = if instrument.is_some() {
            config.precomputed_sparse_scalar_plan.clone()
        } else {
            None
        };
        return Ok(Box::new(PrecomputedTransmissionModel {
            cross_sections: effective_xs,
            density_indices: Arc::new((0..n_params).collect()),
            energies: instrument
                .as_ref()
                .map(|_| Arc::new(config.energies.clone())),
            instrument,
            resolution_plan,
            sparse_cubature_plan,
            sparse_scalar_plan,
            work_layout,
        }));
    }

    let instrument = config
        .resolution
        .clone()
        .map(|r| Arc::new(InstrumentParams { resolution: r }));

    let base_xs = config.precomputed_base_xs.clone();
    let density_ratios = config
        .density_ratios
        .clone()
        .unwrap_or_else(|| vec![1.0; n_density_params]);
    let density_indices = config
        .density_indices
        .clone()
        .unwrap_or_else(|| (0..n_density_params).collect());
    let resolution_plan = if instrument.is_some() {
        config.precomputed_resolution_plan.clone()
    } else {
        None
    };
    let sparse_cubature_plan = if instrument.is_some() {
        config.precomputed_sparse_cubature_plan.clone()
    } else {
        None
    };
    let sparse_scalar_plan = if instrument.is_some() {
        config.precomputed_sparse_scalar_plan.clone()
    } else {
        None
    };
    Ok(Box::new(
        TransmissionFitModel::new(
            config.energies.clone(),
            config.resonance_data.clone(),
            config.temperature_k,
            instrument,
            (density_indices, density_ratios),
            temperature_index,
            base_xs,
        )?
        .with_resolution_plan(resolution_plan)
        .with_sparse_cubature_plan(sparse_cubature_plan)
        .with_sparse_scalar_plan(sparse_scalar_plan),
    ))
}

/// Convert PoissonResult to LmResult with Pearson chi-squared.
fn poisson_to_lm_result(
    model: &dyn FitModel,
    measured_t: &[f64],
    sigma: &[f64],
    pr: &poisson::PoissonResult,
    params: &ParameterSet,
) -> Result<LmResult, PipelineError> {
    let n_free = params.n_free();
    let dof = measured_t.len().saturating_sub(n_free).max(1);
    let y_model = model.evaluate(&pr.params)?;
    let chi_sq: f64 = measured_t
        .iter()
        .zip(y_model.iter())
        .zip(sigma.iter())
        .map(|((obs, mdl), s)| {
            let residual = obs - mdl;
            (residual * residual) / (s * s).max(1e-30)
        })
        .sum();
    Ok(LmResult {
        chi_squared: chi_sq,
        reduced_chi_squared: chi_sq / dof as f64,
        iterations: pr.iterations,
        converged: pr.converged,
        params: pr.params.clone(),
        covariance: pr.covariance.clone(),
        uncertainties: pr.uncertainties.clone(),
    })
}

/// Extract SpectrumFitResult from solver output.
fn extract_result(
    config: &UnifiedFitConfig,
    result: &LmResult,
    n_density_params: usize,
    bg_indices: Option<BackgroundIndices>,
) -> Result<SpectrumFitResult, PipelineError> {
    let densities: Vec<f64> = (0..n_density_params).map(|i| result.params[i]).collect();

    let (anorm, background, back_d, back_f): (f64, [f64; 3], Option<f64>, Option<f64>) =
        if let Some(bi) = bg_indices {
            // `Option<f64>` distinguishes "exponential tail not fit"
            // (`None`) from "fit produced zero" (`Some(0.0)`).
            let bd = bi.back_d.map(|i| result.params[i]);
            let bf = bi.back_f.map(|i| result.params[i]);
            (
                result.params[bi.anorm],
                [
                    result.params[bi.back_a],
                    result.params[bi.back_b],
                    result.params[bi.back_c],
                ],
                bd,
                bf,
            )
        } else {
            (1.0, [0.0, 0.0, 0.0], None, None)
        };

    let (uncertainties, temperature_k, temperature_k_unc) = if result.converged {
        match &result.uncertainties {
            Some(unc_all) => {
                let (temp_k, temp_unc) = if config.fit_temperature {
                    (
                        Some(result.params[n_density_params]),
                        Some(*unc_all.get(n_density_params).unwrap_or(&f64::NAN)),
                    )
                } else {
                    (None, None)
                };
                let unc = unc_all
                    .get(..n_density_params)
                    .map(|s| s.to_vec())
                    .unwrap_or_else(|| vec![f64::NAN; n_density_params]);
                (Some(unc), temp_k, temp_unc)
            }
            None => {
                let temp_k = if config.fit_temperature {
                    Some(result.params[n_density_params])
                } else {
                    None
                };
                (None, temp_k, None)
            }
        }
    } else {
        let temp_k = if config.fit_temperature {
            Some(result.params[n_density_params])
        } else {
            None
        };
        (None, temp_k, None)
    };

    Ok(SpectrumFitResult {
        densities,
        uncertainties,
        reduced_chi_squared: result.reduced_chi_squared,
        converged: result.converged,
        iterations: result.iterations,
        temperature_k,
        temperature_k_unc,
        anorm,
        background,
        back_d,
        back_f,
        t0_us: None,
        l_scale: None,
        deviance_per_dof: None,
    })
}

// ── Research: exact Jacobian/Fisher at arbitrary parameters ──────────────

/// Result of Jacobian/Fisher evaluation at given parameters.
///
/// Produced by [`evaluate_jacobian_and_fisher`], which builds the same model
/// chain as the production fitting pipeline but evaluates at the caller's
/// parameter values instead of optimising.
pub struct ModelJacobianResult {
    /// Analytical Jacobian J (n_data × n_free), row-major.
    pub jacobian: lm::FlatMatrix,
    /// Expected Poisson Fisher F = Jᵀ diag(1/μ) J (n_free × n_free).
    pub fisher: lm::FlatMatrix,
    /// Model prediction μ(E) at the evaluation point.
    pub model_prediction: Vec<f64>,
    /// Names of free parameters, in Jacobian column order.
    pub param_names: Vec<String>,
}

/// Evaluate the exact resolved analytical Jacobian and expected Poisson Fisher
/// at given parameter values, using the same model construction as the
/// production counts-domain fitting pipeline.
///
/// This is a research-oriented function: it builds the full model chain
/// (transmission model → optional background wrappers → counts model),
/// evaluates once at the provided parameters, computes the analytical
/// Jacobian, and assembles the expected Fisher information matrix.
///
/// No optimisation is performed.
///
/// # Arguments
///
/// * `config` — Unified fit configuration (energies, resonance data,
///   resolution, initial_densities used as evaluation densities, etc.)
/// * `flux` — Open-beam counts Φ(E) (length = n_energy)
/// * `background` — Detector background B(E) (length = n_energy, zeros if none)
///
/// Density evaluation values come from `config.initial_densities`.
/// Temperature evaluation value comes from `config.temperature_k`.
/// α₁/α₂ evaluation values come from `config.counts_background` init fields.
///
/// # Errors
/// Returns [`PipelineError::ShapeMismatch`] if `config.precomputed_cross_sections`
/// is set but malformed (empty, wrong row count, wrong row length, or contains a
/// non-finite σ), consistent with `fit_spectrum_typed` and `spatial_map_typed`.
pub fn evaluate_jacobian_and_fisher(
    config: &UnifiedFitConfig,
    flux: &[f64],
    background: &[f64],
) -> Result<ModelJacobianResult, PipelineError> {
    // Reject a malformed caller-supplied precomputed cross-section stack here,
    // before it reaches `build_transmission_model`.  Without this, an empty
    // stack is caught only deep in `PrecomputedTransmissionModel` (as an
    // `InvalidConfig`, not a `ShapeMismatch`) and a wrong-shaped / non-finite
    // stack indexes `xs[0]` / `neg_opt[j]` directly.  This is the third public
    // entry point that consumes `config.precomputed_cross_sections`; it must
    // guard the same way as `fit_spectrum_typed` / `spatial_map_typed` so all
    // three surface the same typed `ShapeMismatch` at the boundary.
    validate_precomputed_cross_sections(config)?;

    let n_density_params = config.n_density_params();

    // ── Build parameter vector — preserves the pre-collapse fixed-flux
    // layout (density → temperature → transmission_background → α₁/α₂)
    // that the research Fisher helper has historically used.  NOT
    // equivalent to the production counts-KL dispatch's parameter
    // layout; kept unchanged for Epic #394 research-script compatibility.
    //     ─────────────────────────────────────────────────────────
    let mut param_vec = build_density_params(config);

    let temperature_index = if config.fit_temperature {
        let idx = param_vec.len();
        param_vec.push(FitParameter {
            name: "temperature_k".into(),
            value: config.temperature_k,
            lower: 1.0,
            upper: 5000.0,
            fixed: false,
        });
        Some(idx)
    } else {
        None
    };

    let kl_bg = if config.transmission_background.is_some() {
        let base = param_vec.len();
        param_vec.push(FitParameter {
            name: "kl_b0".into(),
            value: 0.0,
            lower: 0.0,
            upper: 0.5,
            fixed: false,
        });
        param_vec.push(FitParameter {
            name: "kl_b1".into(),
            value: 0.0,
            lower: 0.0,
            upper: 0.5,
            fixed: false,
        });
        Some((base, base + 1))
    } else {
        None
    };

    let counts_bg = if let Some(bg) = config.counts_background() {
        let alpha1_idx = param_vec.len();
        param_vec.push(if bg.fit_alpha_1 {
            FitParameter {
                name: "alpha_1".into(),
                value: bg.alpha_1_init,
                lower: 0.0,
                upper: 10.0,
                fixed: false,
            }
        } else {
            FitParameter::fixed("alpha_1", bg.alpha_1_init)
        });
        let alpha2_idx = param_vec.len();
        param_vec.push(if bg.fit_alpha_2 {
            FitParameter {
                name: "alpha_2".into(),
                value: bg.alpha_2_init,
                lower: 0.0,
                upper: 10.0,
                fixed: false,
            }
        } else {
            FitParameter::fixed("alpha_2", bg.alpha_2_init)
        });
        Some((alpha1_idx, alpha2_idx))
    } else {
        None
    };

    let params = ParameterSet::new(param_vec);
    let all_vals = params.all_values();
    let free_idx = params.free_indices();
    let n_free = free_idx.len();

    // Collect free parameter names.
    let param_names: Vec<String> = free_idx
        .iter()
        .map(|&i| params.params[i].name.to_string())
        .collect();

    // ── Precompute cross-sections so that analytical Jacobian is available ──
    // For the density-only case (no temperature fitting), the model uses
    // PrecomputedTransmissionModel which requires precomputed XS.
    // For the temperature case, TransmissionFitModel computes base_xs in
    // its constructor.  Either way, precomputing here ensures the analytical
    // Jacobian path is always available.
    // Issue #608: the non-temperature path uses `PrecomputedTransmissionModel`,
    // which must apply resolution on the WORKING grid (auxiliary extended grid
    // under Gaussian resolution) and extract the data points last — the same
    // #608 path as production fitting + spatial mapping, not the old coarse
    // data-grid path.  Derive σ from `resonance_data` (the source of truth) on
    // the working grid here whenever it is not already set up:
    //   * `fit_temperature` → skip: `TransmissionFitModel` builds its own
    //     working-grid base σ internally.
    //   * working-grid σ already attached (a caller did the #608 setup) → skip.
    //   * otherwise (caller passed no σ, OR only a data-grid σ) → (re)build the
    //     data-grid σ = `extract(work σ)` AND the working-grid σ + layout from
    //     `resonance_data`.  A malformed caller-supplied data-grid σ has already
    //     been rejected by the up-front `validate_precomputed_cross_sections`;
    //     here a valid caller σ is superseded by the resonance-data-derived
    //     working-grid σ — the only #608-correct source under Gaussian
    //     resolution.  Without this, a caller that pre-supplied a data-grid σ
    //     plus Gaussian resolution silently got coarse-grid broadening in the
    //     Jacobian/Fisher (R3 finding).
    let config_with_xs;
    let effective_config =
        if config.fit_temperature || config.precomputed_work_cross_sections.is_some() {
            config
        } else {
            let instrument = config
                .resolution
                .clone()
                .map(|r| Arc::new(InstrumentParams { resolution: r }));
            let working = nereids_physics::transmission::broadened_cross_sections_on_working_grid(
                config.energies(),
                &config.resonance_data,
                config.temperature_k,
                instrument.as_deref(),
                None,
            )
            .map_err(PipelineError::Transmission)?;
            config_with_xs = if working.layout.is_identity() {
                // Tabulated / no resolution: the working grid IS the data grid.
                config
                    .clone()
                    .with_precomputed_cross_sections(Arc::new(working.sigma))
            } else {
                // Gaussian aux grid: attach BOTH the extracted data-grid σ (for the
                // surrogate-plan builders + shape validation) and the working-grid σ
                // + layout (AFTER `with_precomputed_cross_sections`, which clears any
                // stale work σ).
                let data_xs: Vec<Vec<f64>> = working
                    .sigma
                    .iter()
                    .map(|s| working.layout.extract(s))
                    .collect();
                config
                    .clone()
                    .with_precomputed_cross_sections(Arc::new(data_xs))
                    .with_precomputed_work_cross_sections(
                        Arc::new(working.sigma),
                        Arc::new(working.layout),
                    )
            };
            &config_with_xs
        };

    // ── Build transmission model (same as production path) ──────────
    let t_model = build_transmission_model(effective_config, n_density_params, temperature_index)?;

    // ── Build counts model chain and evaluate ───────────────────────
    // Use a closure that evaluates and computes Jacobian for any FitModel.
    let evaluate_and_jacobian =
        |model: &dyn FitModel| -> Result<(Vec<f64>, lm::FlatMatrix), PipelineError> {
            let y_model = model.evaluate(&all_vals)?;
            let jac = model
                .analytical_jacobian(&all_vals, &free_idx, &y_model)
                .ok_or_else(|| {
                    PipelineError::InvalidParameter(
                        "analytical Jacobian not available for this model configuration".into(),
                    )
                })?;
            Ok((y_model, jac))
        };

    let (y_model, jac) = if let Some((b0_idx, b1_idx)) = kl_bg {
        let inv_sqrt_e: Vec<f64> = config
            .energies()
            .iter()
            .map(|&e| 1.0 / e.max(1e-10).sqrt())
            .collect();
        let wrapped = poisson::TransmissionKLBackgroundModel {
            inner: &*t_model,
            inv_sqrt_energies: inv_sqrt_e,
            b0_index: b0_idx,
            b1_index: b1_idx,
            n_params: params.params.len(),
        };
        if let Some((a1, a2)) = counts_bg {
            let cm = poisson::CountsBackgroundScaleModel {
                transmission_model: &wrapped,
                flux,
                background,
                alpha1_index: a1,
                alpha2_index: a2,
                n_params: params.params.len(),
            };
            evaluate_and_jacobian(&cm)?
        } else {
            let cm = poisson::CountsModel {
                transmission_model: &wrapped,
                flux,
                background,
                n_params: params.params.len(),
            };
            evaluate_and_jacobian(&cm)?
        }
    } else if let Some((a1, a2)) = counts_bg {
        let cm = poisson::CountsBackgroundScaleModel {
            transmission_model: &*t_model,
            flux,
            background,
            alpha1_index: a1,
            alpha2_index: a2,
            n_params: params.params.len(),
        };
        evaluate_and_jacobian(&cm)?
    } else {
        let cm = poisson::CountsModel {
            transmission_model: &*t_model,
            flux,
            background,
            n_params: params.params.len(),
        };
        evaluate_and_jacobian(&cm)?
    };

    // ── Assemble expected Poisson Fisher: F = Jᵀ diag(1/μ) J ───────
    let mut fisher = lm::FlatMatrix::zeros(n_free, n_free);
    for (i, &mu_i) in y_model.iter().enumerate() {
        let mu_inv = 1.0 / mu_i.max(1e-30);
        for a in 0..n_free {
            let ja = jac.get(i, a);
            for b in 0..=a {
                let jb = jac.get(i, b);
                *fisher.get_mut(a, b) += ja * jb * mu_inv;
                if a != b {
                    *fisher.get_mut(b, a) += ja * jb * mu_inv;
                }
            }
        }
    }

    Ok(ModelJacobianResult {
        jacobian: jac,
        fisher,
        model_prediction: y_model,
        param_names,
    })
}

// ── End Phase 2 ──────────────────────────────────────────────────────────

/// Errors from `FitConfig` construction.
#[derive(Debug, PartialEq)]
pub enum FitConfigError {
    /// Energy grid must be non-empty.
    EmptyEnergies,
    /// Resonance data must be non-empty.
    EmptyResonanceData,
    /// initial_densities length must match resonance_data length.
    DensityCountMismatch { densities: usize, isotopes: usize },
    /// isotope_names length must match resonance_data length.
    NameCountMismatch { names: usize, isotopes: usize },
    /// resonance_data count must match group member count.
    GroupMemberCountMismatch {
        group_name: String,
        rd_count: usize,
        member_count: usize,
    },
    /// ResonanceData isotope doesn't match expected group member.
    GroupMemberIsotopeMismatch {
        group_name: String,
        expected_z: u32,
        expected_a: u32,
        got_z: u32,
        got_a: u32,
    },
    /// Temperature must be finite.
    NonFiniteTemperature(f64),
    /// Temperature must be non-negative.
    NegativeTemperature(f64),
    /// `fit_energy_range` bounds were non-finite or reversed/empty.
    InvalidFitEnergyRange(&'static str),
}

impl fmt::Display for FitConfigError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::EmptyEnergies => write!(f, "energy grid must be non-empty"),
            Self::EmptyResonanceData => write!(f, "resonance_data must be non-empty"),
            Self::DensityCountMismatch {
                densities,
                isotopes,
            } => write!(
                f,
                "initial_densities length ({densities}) must match number of density parameters ({isotopes})"
            ),
            Self::NameCountMismatch { names, isotopes } => write!(
                f,
                "isotope_names length ({names}) must match resonance_data length ({isotopes})"
            ),
            Self::GroupMemberCountMismatch {
                group_name,
                rd_count,
                member_count,
            } => write!(
                f,
                "group '{group_name}': provided {rd_count} ResonanceData but group has {member_count} members"
            ),
            Self::GroupMemberIsotopeMismatch {
                group_name,
                expected_z,
                expected_a,
                got_z,
                got_a,
            } => write!(
                f,
                "group '{group_name}': expected Z={expected_z} A={expected_a} but got Z={got_z} A={got_a}"
            ),
            Self::NonFiniteTemperature(v) => {
                write!(f, "temperature must be finite, got {v}")
            }
            Self::NegativeTemperature(v) => {
                write!(f, "temperature must be non-negative, got {v}")
            }
            Self::InvalidFitEnergyRange(msg) => {
                write!(f, "invalid fit_energy_range: {msg}")
            }
        }
    }
}

impl std::error::Error for FitConfigError {}

/// Result of fitting a single spectrum.
#[derive(Debug, Clone)]
pub struct SpectrumFitResult {
    /// Fitted areal densities (atoms/barn), one per isotope.
    pub densities: Vec<f64>,
    /// Uncertainty on each density.
    ///
    /// `None` when covariance computation was skipped.
    pub uncertainties: Option<Vec<f64>>,
    /// Reduced chi-squared of the fit.
    pub reduced_chi_squared: f64,
    /// Whether the fit converged.
    pub converged: bool,
    /// Number of iterations.
    pub iterations: usize,
    /// Fitted temperature in Kelvin (only when `fit_temperature` is true).
    pub temperature_k: Option<f64>,
    /// 1-sigma uncertainty on the fitted temperature (from covariance matrix).
    pub temperature_k_unc: Option<f64>,
    /// Fitted normalization scale (SAMMY `Anorm`).  When background
    /// fitting is disabled, the pipeline emits `1.0` (memo 35 §P1
    /// convention — `λ̂` absorbs the scale).
    pub anorm: f64,
    /// Fitted background polynomial coefficients `[BackA, BackB, BackC]`
    /// for the SAMMY-style 6-term background:
    ///
    /// ```text
    /// bg(E) = BackA + BackB / √E + BackC · √E + BackD · exp(-BackF / √E)
    /// ```
    ///
    /// Both LM-transmission and counts-KL solver paths use the same
    /// SAMMY semantics here — the legacy alpha-fitting `[b0, b1,
    /// alpha_2]` layout was removed when `fit_counts_poisson` was
    /// retired in PR #450.  When background fitting is disabled, the
    /// pipeline emits `[0.0, 0.0, 0.0]`.
    pub background: [f64; 3],
    /// Fitted exponential background amplitude (SAMMY BackD).
    /// `None` when the exponential tail is not fitted; `Some(value)`
    /// when the LM transmission background was active with
    /// `fit_back_d=true`.  Mirrors [`Self::t0_us`] /
    /// [`Self::l_scale`] semantics so the GUI overlay can
    /// distinguish "unfit" from "fitted to zero" without an ambiguous
    /// `0.0` sentinel.
    pub back_d: Option<f64>,
    /// Fitted exponential background decay constant (SAMMY BackF).
    /// `None` when the exponential tail is not fitted; `Some(value)`
    /// when the LM transmission background was active with
    /// `fit_back_f=true`.  Paired with [`Self::back_d`] — SAMMY
    /// requires both flags toggled together (see
    /// `validate_transmission_background`).
    pub back_f: Option<f64>,
    /// Fitted TOF offset in microseconds (SAMMY TZERO t₀).
    /// `None` when energy-scale fitting is not enabled.
    pub t0_us: Option<f64>,
    /// Fitted flight-path scale factor (SAMMY TZERO L₀, dimensionless).
    /// `None` when energy-scale fitting is not enabled.
    pub l_scale: Option<f64>,
    /// Conditional binomial deviance divided by `(n − k)`
    /// (memo 35 §P1.2 — primary GOF for the counts-KL dispatch, i.e.
    /// `SolverConfig::PoissonKL` on `InputData::Counts` or
    /// `InputData::CountsWithNuisance`).
    ///
    /// `Some(D/dof)` when the counts-KL (joint-Poisson) path was used;
    /// `None` for the LM path and for transmission + PoissonKL (those
    /// populate `reduced_chi_squared` with Pearson χ² / (n−k) instead).
    pub deviance_per_dof: Option<f64>,
}

#[cfg(test)]
mod tests {
    use super::*;
    use nereids_endf::resonance::test_support::{
        hf178_mlbw_two_resonances, synthetic_single_resonance, u238_single_resonance,
    };
    use nereids_fitting::lm::FitModel;
    use nereids_physics::transmission as phys_transmission;

    /// Issue #608: the dip detector finds the resonance signatures (local
    /// transmission minima) the energy-scale peak-match seed needs.
    #[test]
    fn detect_transmission_dips_finds_clear_dips() {
        let energies: Vec<f64> = (0..120).map(|i| 1.0 + (i as f64) * 0.5).collect();
        let mut t = vec![1.0_f64; 120];
        // Two clear, multi-bin dips so 3-point smoothing keeps them as minima.
        for v in t.iter_mut().take(33).skip(28) {
            *v = 0.4;
        }
        for v in t.iter_mut().take(83).skip(78) {
            *v = 0.6;
        }
        let dips = detect_transmission_dips(&t, &energies);
        assert_eq!(dips.len(), 2, "expected 2 dips, got {dips:?}");
        let mut de: Vec<f64> = dips.iter().map(|&(e, _)| e).collect();
        de.sort_by(f64::total_cmp);
        // Centres near energies[30] (16.0 eV) and energies[80] (41.0 eV).
        assert!(
            (de[0] - energies[30]).abs() <= 1.0,
            "first dip at {} eV",
            de[0]
        );
        assert!(
            (de[1] - energies[80]).abs() <= 1.0,
            "second dip at {} eV",
            de[1]
        );
    }

    /// Issue #608: the energy-scale peak-match seed recovers the identity
    /// calibration (t0≈0, L_scale≈1) when the transmission dips sit at the known
    /// resonance energies — exercising the full seed path (dip detection,
    /// nearest-resonance matching, linear TOF fit).
    #[test]
    fn peak_match_energy_scale_seed_identity() {
        let data = hf178_mlbw_two_resonances(); // s-waves at 7.8 and 16.9 eV
        let energies: Vec<f64> = (0..400).map(|i| 4.0 + (i as f64) * 0.05).collect();
        let (t_obs, _sigma) = synthetic_transmission(&data, 0.1, &energies);
        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![data],
            vec!["Hf-178".into()],
            293.6,
            None,
            vec![0.1],
        )
        .unwrap()
        .with_energy_scale(0.0, 1.0, 25.0);
        let (t0, l_scale) = peak_match_energy_scale_seed(
            &t_obs,
            config.energies(),
            &config,
            25.0,
            (-10.0, 10.0),
            (0.99, 1.01),
        )
        .expect("seed should be Some with 2 resonances and detectable dips");
        // Dips at the un-shifted resonance energies ⇒ identity calibration.
        assert!(
            t0.abs() < 0.5,
            "t0 should be ≈0 for un-shifted data, got {t0}"
        );
        assert!(
            (l_scale - 1.0).abs() < 5e-3,
            "L_scale should be ≈1 for un-shifted data, got {l_scale}"
        );
    }

    /// Issue #608 (R4): the peak-match seed must recover a NON-identity
    /// calibration, not just confirm identity.  Generate measured data with a
    /// known injected `(t0, L_scale)` via the `EnergyScaleTransmissionModel`
    /// (which shifts the resonance positions), then assert the seed recovers it.
    #[test]
    fn peak_match_energy_scale_seed_recovers_nonidentity() {
        let data = hf178_mlbw_two_resonances(); // s-waves at 7.8 and 16.9 eV
        // Finer grid than the identity test so the discrete dip positions pin
        // the slope (L_scale) tightly enough to distinguish it from 1.0.
        let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
        let density = 0.1_f64;
        let flight_path = 25.0_f64;
        let (t0_true, l_scale_true) = (2.0_f64, 1.006_f64);
        // EnergyScale model (no resolution: dip POSITIONS, not shapes, drive the
        // seed) evaluated at the injected calibration ⇒ shifted measured data.
        let model = EnergyScaleTransmissionModel::new(
            Arc::new(vec![data.clone()]),
            Arc::new(vec![0]),
            Arc::new(vec![1.0]),
            293.6,
            energies.clone(),
            flight_path,
            1, // t0 index
            2, // l_scale index
            None,
        );
        let t_obs = model.evaluate(&[density, t0_true, l_scale_true]).unwrap();
        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![data],
            vec!["Hf-178".into()],
            293.6,
            None,
            vec![density],
        )
        .unwrap()
        .with_energy_scale(0.0, 1.0, flight_path);
        let (t0, l_scale) = peak_match_energy_scale_seed(
            &t_obs,
            config.energies(),
            &config,
            flight_path,
            (-10.0, 10.0),
            (0.99, 1.01),
        )
        .expect("seed should be Some for a clean shifted two-resonance spectrum");
        // The seed is a coarse peak-match (dip energies quantized to the grid),
        // so tolerances are looser than a converged LM — it just needs to land
        // in the global-min basin, clearly distinct from the cold start (0, 1).
        assert!(
            (t0 - t0_true).abs() < 0.6,
            "seed should recover t0 ≈ {t0_true}, got {t0}"
        );
        assert!(
            (l_scale - l_scale_true).abs() < 4e-3,
            "seed should recover L_scale ≈ {l_scale_true}, got {l_scale}"
        );
        // Sanity: the recovered seed is strictly closer to truth than cold start.
        assert!(
            (t0 - t0_true).abs() < (0.0 - t0_true).abs()
                && (l_scale - l_scale_true).abs() < (1.0 - l_scale_true).abs(),
            "seed must improve on the cold start"
        );
    }

    /// Issue #608 (R4): a heavily-absorbing but FEATURELESS spectrum (no
    /// resonance dips) yields fewer than two detectable dips, so the seed must
    /// return `None` and the caller keeps the cold start rather than fitting a
    /// spurious calibration.
    #[test]
    fn peak_match_energy_scale_seed_none_on_featureless_spectrum() {
        let data = hf178_mlbw_two_resonances(); // 2 resonances in range (gate passes)
        let energies: Vec<f64> = (0..400).map(|i| 4.0 + (i as f64) * 0.05).collect();
        // Strong smooth 1/√E absorption, monotonic ⇒ NO local minima ⇒ < 2 dips.
        let t_obs: Vec<f64> = energies.iter().map(|&e| (-5.0 / e.sqrt()).exp()).collect();
        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![data],
            vec!["Hf-178".into()],
            293.6,
            None,
            vec![0.1],
        )
        .unwrap()
        .with_energy_scale(0.0, 1.0, 25.0);
        assert!(
            peak_match_energy_scale_seed(
                &t_obs,
                config.energies(),
                &config,
                25.0,
                (-10.0, 10.0),
                (0.99, 1.01),
            )
            .is_none(),
            "a featureless heavily-absorbing spectrum has no resonance dips ⇒ \
             seed must return None (cold-start fallback)"
        );
    }

    /// Issue #608 (Copilot PR #609): duplicate resonance center energies — two
    /// isotopes with a resonance at the SAME energy in a grouped fit — must not
    /// collapse the seed's `match_tol` to 0.  Pre-fix `resonance_center_energies`
    /// returned `[7.8, 7.8, 16.9]` ⇒ minimum spacing 0 ⇒ `match_tol` 0 ⇒ every
    /// dip rejected ⇒ seed silently `None` (cold start).  With the dedup + grid
    /// floor the seed recovers the (identity) calibration.
    #[test]
    fn peak_match_energy_scale_seed_handles_duplicate_resonance_energies() {
        use nereids_physics::transmission::{SampleParams, forward_model};

        // Two isotopes share an EXACT resonance energy (7.8); a third is distinct.
        let iso_a = synthetic_single_resonance(72, 178, 176.0, 7.8);
        let iso_b = synthetic_single_resonance(74, 184, 182.0, 7.8); // duplicate energy
        let iso_c = synthetic_single_resonance(40, 90, 89.0, 16.9); // distinct
        let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
        let density = 0.05_f64;
        let sample = SampleParams::new(
            293.6,
            vec![
                (iso_a.clone(), density),
                (iso_b.clone(), density),
                (iso_c.clone(), density),
            ],
        )
        .unwrap();
        let t_obs = forward_model(&energies, &sample, None).unwrap();
        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![iso_a, iso_b, iso_c],
            vec!["A".into(), "B".into(), "C".into()],
            293.6,
            None,
            vec![density, density, density],
        )
        .unwrap()
        .with_energy_scale(0.0, 1.0, 25.0);
        let (t0, l_scale) = peak_match_energy_scale_seed(
            &t_obs,
            config.energies(),
            &config,
            25.0,
            (-10.0, 10.0),
            (0.99, 1.01),
        )
        .expect(
            "duplicate resonance energies must not collapse match_tol to 0 — \
             seed should be Some (was silently None pre-fix)",
        );
        assert!(t0.abs() < 0.5, "identity calibration: t0 ≈ 0, got {t0}");
        assert!(
            (l_scale - 1.0).abs() < 5e-3,
            "identity calibration: L_scale ≈ 1, got {l_scale}"
        );
    }

    // ── Phase 0: InputData + SolverConfig + CountsBackgroundConfig tests ──

    #[test]
    fn test_input_data_transmission_n_energies() {
        let data = InputData::Transmission {
            transmission: vec![0.9, 0.8, 0.7],
            uncertainty: vec![0.01, 0.01, 0.01],
        };
        assert_eq!(data.n_energies(), 3);
        assert!(!data.is_counts());
    }

    #[test]
    fn test_input_data_counts_n_energies() {
        let data = InputData::Counts {
            sample_counts: vec![10.0, 20.0, 30.0, 40.0],
            open_beam_counts: vec![100.0, 100.0, 100.0, 100.0],
        };
        assert_eq!(data.n_energies(), 4);
        assert!(data.is_counts());
    }

    #[test]
    fn test_input_data_counts_with_nuisance() {
        let data = InputData::CountsWithNuisance {
            sample_counts: vec![5.0, 6.0],
            flux: vec![100.0, 100.0],
            background: vec![0.5, 0.5],
        };
        assert_eq!(data.n_energies(), 2);
        assert!(data.is_counts());
    }

    #[test]
    fn test_solver_config_default_is_auto() {
        let cfg = SolverConfig::default();
        assert!(matches!(cfg, SolverConfig::Auto));
    }

    #[test]
    fn test_counts_background_config_default() {
        let cfg = CountsBackgroundConfig::default();
        assert!((cfg.alpha_1_init - 1.0).abs() < f64::EPSILON);
        assert!((cfg.alpha_2_init - 1.0).abs() < f64::EPSILON);
        assert!(!cfg.fit_alpha_1);
        assert!(!cfg.fit_alpha_2);
    }

    // ── Phase 2: fit_spectrum_typed tests ──

    /// Helper: build synthetic transmission data from known density.
    fn synthetic_transmission(
        data: &ResonanceData,
        true_density: f64,
        energies: &[f64],
    ) -> (Vec<f64>, Vec<f64>) {
        let model = PrecomputedTransmissionModel {
            cross_sections: Arc::new(vec![
                phys_transmission::broadened_cross_sections(
                    energies,
                    std::slice::from_ref(data),
                    0.0,
                    None,
                    None,
                )
                .unwrap()
                .into_iter()
                .next()
                .unwrap(),
            ]),
            density_indices: Arc::new(vec![0]),
            energies: None,
            instrument: None,
            resolution_plan: None,
            sparse_cubature_plan: None,
            sparse_scalar_plan: None,
            work_layout: None,
        };
        let t = model.evaluate(&[true_density]).unwrap();
        let sigma: Vec<f64> = t.iter().map(|&v| 0.01 * v.max(0.01)).collect();
        (t, sigma)
    }

    /// Helper: build synthetic counts from known density.
    fn synthetic_counts(
        data: &ResonanceData,
        true_density: f64,
        energies: &[f64],
        i0: f64,
    ) -> (Vec<f64>, Vec<f64>) {
        let (t, _) = synthetic_transmission(data, true_density, energies);
        let open_beam: Vec<f64> = vec![i0; energies.len()];
        let sample: Vec<f64> = t.iter().map(|&v| (v * i0).round().max(0.0)).collect();
        (sample, open_beam)
    }

    #[test]
    fn test_typed_transmission_lm_recovers_density() {
        let data = u238_single_resonance();
        let true_density = 0.002;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission(&data, true_density, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged, "LM should converge");
        let fitted = result.densities[0];
        assert!(
            (fitted - true_density).abs() / true_density < 0.05,
            "density: fitted={fitted}, true={true_density}"
        );
    }

    /// Helper: a valid single-isotope transmission config + input on a short
    /// grid, for precomputed-cross-section shape-validation tests.
    fn precomputed_xs_fixture() -> (UnifiedFitConfig, InputData) {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let (t, sigma) = synthetic_transmission(&data, 0.001, &energies);
        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };
        (config, input)
    }

    #[test]
    fn test_fit_rejects_empty_precomputed_cross_sections() {
        // An empty XS stack used to panic on `xs[0].len()` deep in the
        // forward-model builder; it must now be a typed up-front error.
        let (config, input) = precomputed_xs_fixture();
        let config = config.with_precomputed_cross_sections(Arc::new(Vec::new()));
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            matches!(err, PipelineError::ShapeMismatch(_)),
            "expected ShapeMismatch, got {err:?}"
        );
        assert!(err.to_string().contains("must not be empty"));
    }

    #[test]
    fn test_fit_rejects_wrong_row_count_precomputed_cross_sections() {
        // 1 isotope → exactly 1 σ row expected; supply 2.
        let (config, input) = precomputed_xs_fixture();
        let n_e = config.energies().len();
        let bad = vec![vec![1.0; n_e], vec![1.0; n_e]];
        let config = config.with_precomputed_cross_sections(Arc::new(bad));
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            matches!(err, PipelineError::ShapeMismatch(_)),
            "expected ShapeMismatch, got {err:?}"
        );
        assert!(err.to_string().contains("rows"));
    }

    #[test]
    fn test_fit_rejects_wrong_energy_length_precomputed_cross_sections() {
        // Correct row count (1) but the row is the wrong length.
        let (config, input) = precomputed_xs_fixture();
        let n_e = config.energies().len();
        let bad = vec![vec![1.0; n_e + 3]];
        let config = config.with_precomputed_cross_sections(Arc::new(bad));
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            matches!(err, PipelineError::ShapeMismatch(_)),
            "expected ShapeMismatch, got {err:?}"
        );
        assert!(err.to_string().contains("config.energies"));
    }

    #[test]
    fn test_fit_accepts_correct_precomputed_cross_sections() {
        // A correctly-shaped precomputed stack must still fit (no false
        // rejection): 1 row of length n_e.
        let (config, input) = precomputed_xs_fixture();
        let n_e = config.energies().len();
        let xs = phys_transmission::broadened_cross_sections(
            config.energies(),
            config.resonance_data(),
            0.0,
            None,
            None,
        )
        .unwrap();
        assert_eq!(xs.len(), 1);
        assert_eq!(xs[0].len(), n_e);
        let config = config.with_precomputed_cross_sections(Arc::new(xs));
        // Must not error on shape; the fit itself may or may not converge,
        // but the call must reach the solver rather than fail validation.
        let result = fit_spectrum_typed(&input, &config);
        assert!(
            result.is_ok(),
            "correctly-shaped precomputed XS must pass validation, got {result:?}"
        );
    }

    #[test]
    fn test_fit_rejects_non_finite_precomputed_cross_sections() {
        // A correctly-shaped row whose σ values are NaN passes the shape
        // checks but poisons the forward model (NaN residual swallowed as a
        // failed pixel).  It must be rejected at the boundary.  `NaN < x` is
        // `false`, so a bare order comparison would let it through — the
        // `is_finite()` half of the guard is what catches it.
        let (config, input) = precomputed_xs_fixture();
        let n_e = config.energies().len();
        let bad = vec![vec![f64::NAN; n_e]];
        let config = config.with_precomputed_cross_sections(Arc::new(bad));
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            matches!(err, PipelineError::ShapeMismatch(_)),
            "expected ShapeMismatch, got {err:?}"
        );
        assert!(
            err.to_string().contains("non-finite"),
            "error should mention non-finite σ, got: {err}"
        );

        // ±∞ is equally rejected.
        let (config, input) = precomputed_xs_fixture();
        let mut row = vec![1.0; n_e];
        row[n_e / 2] = f64::INFINITY;
        let config = config.with_precomputed_cross_sections(Arc::new(vec![row]));
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            matches!(err, PipelineError::ShapeMismatch(_)),
            "expected ShapeMismatch for +inf σ, got {err:?}"
        );
    }

    #[test]
    fn test_evaluate_jacobian_rejects_malformed_precomputed_cross_sections() {
        // `evaluate_jacobian_and_fisher` is a third public entry point that
        // consumes `config.precomputed_cross_sections`.  An empty stack used
        // to reach `build_transmission_model` and panic on `xs[0].len()`; it
        // must now fail the shared up-front validator instead.
        let (config, _input) = precomputed_xs_fixture();
        let n_e = config.energies().len();
        let flux = vec![1.0; n_e];
        let background = vec![0.0; n_e];

        // `ModelJacobianResult` (the Ok type) is not `Debug`, so match on the
        // result rather than calling `unwrap_err()`.
        let empty = config
            .clone()
            .with_precomputed_cross_sections(Arc::new(Vec::new()));
        match evaluate_jacobian_and_fisher(&empty, &flux, &background) {
            Err(PipelineError::ShapeMismatch(msg)) => {
                assert!(msg.contains("must not be empty"), "got: {msg}");
            }
            Err(other) => panic!("expected ShapeMismatch for empty XS, got {other:?}"),
            Ok(_) => panic!("empty precomputed XS must be rejected"),
        }

        // Non-finite σ is rejected on this path too.
        let nan = config.with_precomputed_cross_sections(Arc::new(vec![vec![f64::NAN; n_e]]));
        match evaluate_jacobian_and_fisher(&nan, &flux, &background) {
            Err(PipelineError::ShapeMismatch(msg)) => {
                assert!(msg.contains("non-finite"), "got: {msg}");
            }
            Err(other) => panic!("expected ShapeMismatch for NaN σ, got {other:?}"),
            Ok(_) => panic!("non-finite precomputed σ must be rejected"),
        }
    }

    #[test]
    fn test_extract_result_drops_uncertainties_when_unconverged() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            293.6,
            None,
            vec![0.001],
        )
        .unwrap();

        let result = LmResult {
            chi_squared: 1.0,
            reduced_chi_squared: 1.0,
            iterations: 5,
            converged: false,
            params: vec![0.001],
            covariance: Some(lm::FlatMatrix::zeros(1, 1)),
            uncertainties: Some(vec![0.123]),
        };

        let extracted = extract_result(&config, &result, 1, None).unwrap();
        assert!(!extracted.converged);
        assert!(
            extracted.uncertainties.is_none(),
            "pipeline must not surface uncertainties from an unconverged fit"
        );
        assert!(extracted.temperature_k_unc.is_none());
    }

    #[test]
    fn test_typed_counts_kl_recovers_density() {
        let data = u238_single_resonance();
        let true_density = 0.002;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged, "KL on counts should converge");
        let fitted = result.densities[0];
        assert!(
            (fitted - true_density).abs() / true_density < 0.10,
            "density: fitted={fitted}, true={true_density}"
        );
    }

    #[test]
    fn test_typed_counts_kl_low_counts_recovers_density() {
        // I0=10 counts per bin — the regime where KL excels and LM fails
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 10.0);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged, "KL on low counts should converge");
        let fitted = result.densities[0];
        // Wider tolerance for low counts
        assert!(
            (fitted - true_density).abs() / true_density < 0.30,
            "density: fitted={fitted}, true={true_density}"
        );
    }

    #[test]
    fn test_typed_transmission_kl_recovers_density() {
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission(&data, true_density, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged, "KL on transmission should converge");
        let fitted = result.densities[0];
        assert!(
            (fitted - true_density).abs() / true_density < 0.05,
            "density: fitted={fitted}, true={true_density}"
        );
    }

    #[test]
    fn test_typed_counts_lm_auto_converts() {
        // Counts + LM should auto-convert to transmission and fit
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(
            result.converged,
            "LM on auto-converted counts should converge"
        );
        let fitted = result.densities[0];
        assert!(
            (fitted - true_density).abs() / true_density < 0.10,
            "density: fitted={fitted}, true={true_density}"
        );
    }

    #[test]
    fn test_typed_auto_solver_selects_kl_for_counts() {
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);

        // Auto solver (default)
        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap(); // SolverConfig::Auto by default

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(
            result.converged,
            "Auto solver on counts should use KL and converge"
        );
    }

    #[test]
    fn test_typed_transmission_with_background_lm() {
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let true_anorm = 0.95;
        let true_back_a = 0.02;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();

        // Generate synthetic data with background
        let (t_pure, _) = synthetic_transmission(&data, true_density, &energies);
        let t_bg: Vec<f64> = t_pure
            .iter()
            .map(|&v| true_anorm * v + true_back_a)
            .collect();
        let sigma: Vec<f64> = t_bg.iter().map(|&v| 0.01 * v.max(0.01)).collect();

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig {
            max_iter: 500,
            ..LmConfig::default()
        }))
        .with_transmission_background(BackgroundConfig::default());

        let input = InputData::Transmission {
            transmission: t_bg,
            uncertainty: sigma,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(
            result.converged,
            "LM+BG should converge on noiseless synthetic data (chi2r={}, iter={})",
            result.reduced_chi_squared, result.iterations
        );
        assert!(
            (result.densities[0] - true_density).abs() / true_density < 0.05,
            "density: fitted={}, true={true_density}",
            result.densities[0]
        );
        assert!(
            (result.anorm - true_anorm).abs() / true_anorm < 0.05,
            "anorm: fitted={}, true={true_anorm}",
            result.anorm
        );
    }

    #[test]
    fn test_typed_counts_with_nuisance_rejects_lm() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..10).map(|i| 1.0 + (i as f64) * 0.5).collect();

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));

        let input = InputData::CountsWithNuisance {
            sample_counts: vec![10.0; 10],
            flux: vec![100.0; 10],
            background: vec![0.0; 10],
        };

        let result = fit_spectrum_typed(&input, &config);
        assert!(result.is_err(), "CountsWithNuisance + LM should error");
    }

    /// Helper: build synthetic transmission at a given temperature.
    fn synthetic_transmission_at_temp(
        data: &ResonanceData,
        true_density: f64,
        temperature_k: f64,
        energies: &[f64],
    ) -> (Vec<f64>, Vec<f64>) {
        let xs = phys_transmission::broadened_cross_sections(
            energies,
            std::slice::from_ref(data),
            temperature_k,
            None,
            None,
        )
        .unwrap();
        let model = PrecomputedTransmissionModel {
            cross_sections: Arc::new(xs),
            density_indices: Arc::new(vec![0]),
            energies: None,
            instrument: None,
            resolution_plan: None,
            sparse_cubature_plan: None,
            sparse_scalar_plan: None,
            work_layout: None,
        };
        let t = model.evaluate(&[true_density]).unwrap();
        let sigma: Vec<f64> = t.iter().map(|&v| 0.01 * v.max(0.01)).collect();
        (t, sigma)
    }

    #[test]
    fn test_typed_poisson_kl_with_temperature() {
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let true_temp = 350.0;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            300.0, // initial guess (off by 50 K)
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig {
            max_iter: 500,
            ..PoissonConfig::default()
        }))
        .with_fit_temperature(true);

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();

        // Check density recovery (within 1%)
        let fitted_density = result.densities[0];
        assert!(
            (fitted_density - true_density).abs() / true_density < 0.01,
            "density: fitted={fitted_density}, true={true_density}, ratio={}",
            (fitted_density - true_density).abs() / true_density,
        );

        // Check temperature recovery (within 1 K)
        let fitted_temp = result
            .temperature_k
            .expect("temperature_k should be Some when fit_temperature=true");
        assert!(
            (fitted_temp - true_temp).abs() < 1.0,
            "temperature: fitted={fitted_temp}, true={true_temp}, delta={}",
            (fitted_temp - true_temp).abs(),
        );
    }

    #[test]
    fn test_typed_poisson_kl_with_temperature_and_background() {
        let data = u238_single_resonance();
        let true_density = 0.0005;
        let true_temp = 350.0;
        let true_b0 = 0.012;
        let true_b1 = 0.008;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission_at_temp(&data, true_density, true_temp, &energies);
        let measured_t: Vec<f64> = t
            .iter()
            .zip(energies.iter())
            .map(|(&ti, &e)| ti + true_b0 + true_b1 / e.sqrt())
            .collect();

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            300.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig {
            max_iter: 120,
            gauss_newton_lambda: 1e-4,
            ..PoissonConfig::default()
        }))
        .with_fit_temperature(true)
        .with_transmission_background(BackgroundConfig::default());

        let input = InputData::Transmission {
            transmission: measured_t,
            uncertainty: sigma,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();

        assert!(result.converged, "fit did not converge: {result:?}");
        assert!(
            result.iterations <= 80,
            "expected KL background+temperature fit to converge well before max_iter; got {}",
            result.iterations,
        );

        let fitted_density = result.densities[0];
        assert!(
            (fitted_density - true_density).abs() / true_density < 0.02,
            "density: fitted={fitted_density}, true={true_density}, ratio={}",
            (fitted_density - true_density).abs() / true_density,
        );

        let fitted_temp = result
            .temperature_k
            .expect("temperature_k should be Some when fit_temperature=true");
        assert!(
            (fitted_temp - true_temp).abs() < 3.0,
            "temperature: fitted={fitted_temp}, true={true_temp}, delta={}",
            (fitted_temp - true_temp).abs(),
        );

        assert!(
            (result.background[0] - true_b0).abs() < 5e-3,
            "background b0: fitted={}, true={}",
            result.background[0],
            true_b0,
        );
        assert!(
            (result.background[1] - true_b1).abs() < 5e-3,
            "background b1: fitted={}, true={}",
            result.background[1],
            true_b1,
        );
    }

    /// Round-trip test: create a group of 2 isotopes with known ratios,
    /// generate synthetic transmission, fit with group constraints,
    /// verify the fitted group density matches the true value.
    #[test]
    fn test_grouped_fit_spectrum_round_trip() {
        use nereids_core::types::IsotopeGroup;

        // Two synthetic isotopes with resonances at different energies
        let rd1 = synthetic_single_resonance(92, 235, 233.025, 5.0);
        let rd2 = synthetic_single_resonance(92, 238, 236.006, 7.0);

        // Group with 60/40 ratio
        let iso1 = nereids_core::types::Isotope::new(92, 235).unwrap();
        let iso2 = nereids_core::types::Isotope::new(92, 238).unwrap();
        let group =
            IsotopeGroup::custom("U (60/40)".into(), vec![(iso1, 0.6), (iso2, 0.4)]).unwrap();

        let energies: Vec<f64> = (0..301).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let true_density = 0.0005;

        // Generate synthetic transmission using effective densities
        let sample = nereids_physics::transmission::SampleParams::new(
            0.0,
            vec![
                (rd1.clone(), true_density * 0.6),
                (rd2.clone(), true_density * 0.4),
            ],
        )
        .unwrap();
        let transmission =
            nereids_physics::transmission::forward_model(&energies, &sample, None).unwrap();
        let uncertainty: Vec<f64> = transmission.iter().map(|&t| 0.01 * t.max(0.01)).collect();

        // Build config with group
        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![rd1.clone()],
            vec!["placeholder".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_groups(&[(&group, &[rd1, rd2])], vec![0.001])
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));

        let input = InputData::Transmission {
            transmission,
            uncertainty,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();

        // Should recover true density within 1%
        assert_eq!(result.densities.len(), 1, "should have 1 group density");
        let fitted = result.densities[0];
        let rel_error = (fitted - true_density).abs() / true_density;
        assert!(
            rel_error < 0.01,
            "group density: fitted={fitted}, true={true_density}, rel_error={rel_error}"
        );
        assert!(result.converged, "fit should converge");
    }

    #[test]
    fn test_grouped_poisson_kl_with_temperature_and_background_noiseless() {
        use nereids_core::types::IsotopeGroup;

        let rd1 = synthetic_single_resonance(72, 176, 8.5, 5.0);
        let rd2 = synthetic_single_resonance(72, 178, 17.0, 7.5);
        let rd3 = synthetic_single_resonance(72, 180, 29.0, 6.0);

        let hf176 = nereids_core::types::Isotope::new(72, 176).unwrap();
        let hf178 = nereids_core::types::Isotope::new(72, 178).unwrap();
        let hf180 = nereids_core::types::Isotope::new(72, 180).unwrap();
        let group = IsotopeGroup::custom(
            "Hf-like (3 member)".into(),
            vec![(hf176, 0.2), (hf178, 0.5), (hf180, 0.3)],
        )
        .unwrap();

        let energies: Vec<f64> = (0..300).map(|i| 1.0 + (49.0 * i as f64) / 299.0).collect();
        let true_density = 0.001;
        let true_temp = 400.0;
        let true_b0 = 0.012;
        let true_b1 = 0.008;

        let sample = nereids_physics::transmission::SampleParams::new(
            true_temp,
            vec![
                (rd1.clone(), true_density * 0.2),
                (rd2.clone(), true_density * 0.5),
                (rd3.clone(), true_density * 0.3),
            ],
        )
        .unwrap();
        let pure_t =
            nereids_physics::transmission::forward_model(&energies, &sample, None).unwrap();
        let measured_t: Vec<f64> = pure_t
            .iter()
            .zip(energies.iter())
            .map(|(&t, &e)| t + true_b0 + true_b1 / e.sqrt())
            .collect();
        let sigma = vec![0.001; energies.len()];

        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![rd1.clone()],
            vec!["placeholder".into()],
            293.6,
            None,
            vec![0.0008],
        )
        .unwrap()
        .with_groups(&[(&group, &[rd1, rd2, rd3])], vec![0.0008])
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig {
            max_iter: 200,
            gauss_newton_lambda: 1e-4,
            ..PoissonConfig::default()
        }))
        .with_fit_temperature(true)
        .with_transmission_background(BackgroundConfig::default());

        let input = InputData::Transmission {
            transmission: measured_t,
            uncertainty: sigma,
        };

        let result = fit_spectrum_typed(&input, &config).unwrap();

        assert!(result.converged, "fit did not converge: {result:?}");
        // Tolerance: 1% — the NormalizedTransmissionModel (4 background params)
        // has slightly different convergence than the old 2-param KL model.
        assert!(
            (result.densities[0] - true_density).abs() / true_density < 0.01,
            "density: fitted={}, true={true_density}",
            result.densities[0]
        );
        let fitted_temp = result
            .temperature_k
            .expect("temperature_k should be Some when fit_temperature=true");
        assert!(
            (fitted_temp - true_temp).abs() < 8.0,
            "temperature: fitted={fitted_temp}, true={true_temp}",
        );
        // Background: NormalizedTransmissionModel distributes the additive
        // background across Anorm + BackA/B/C.  Check that the total background
        // contribution is reasonable, not individual parameters.
        let e_mid: f64 = 10.0;
        let bg_total = (result.anorm - 1.0)
            + result.background[0]
            + result.background[1] / e_mid.sqrt()
            + result.background[2] * e_mid.sqrt();
        let true_bg_mid = true_b0 + true_b1 / e_mid.sqrt();
        assert!(
            (bg_total - true_bg_mid).abs() < 0.02,
            "total bg at E={e_mid}: fitted={bg_total:.6}, true={true_bg_mid:.6}",
        );
    }

    // ── Phase 2: KL fitting uncertainty tests ──────────────────────────────

    #[test]
    fn test_kl_counts_returns_density_uncertainty() {
        let data = u238_single_resonance();
        let true_density = 0.002;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, true_density, &energies, 1000.0);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged);
        let unc = result
            .uncertainties
            .as_ref()
            .expect("KL 1D fit should return density uncertainties");
        assert_eq!(unc.len(), 1);
        assert!(
            unc[0].is_finite() && unc[0] > 0.0,
            "density unc = {}",
            unc[0]
        );
        assert!(
            unc[0] < result.densities[0],
            "unc ({}) should be < density ({}) for high-count data",
            unc[0],
            result.densities[0]
        );
    }

    #[test]
    fn test_kl_counts_returns_temperature_uncertainty() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            350.0,
            None,
            vec![0.0005],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_fit_temperature(true);

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged);
        let unc = result
            .uncertainties
            .as_ref()
            .expect("KL+temp fit should return density uncertainties");
        assert!(
            unc[0].is_finite() && unc[0] > 0.0,
            "density unc = {}",
            unc[0]
        );
        let t_unc = result
            .temperature_k_unc
            .expect("KL+temp fit should return temperature uncertainty");
        assert!(
            t_unc.is_finite() && t_unc > 0.0,
            "temperature unc = {t_unc}"
        );
    }

    #[test]
    fn test_kl_counts_with_background_returns_uncertainty() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..201).map(|i| 4.0 + (i as f64) * 0.05).collect();
        let (sample, open_beam) = synthetic_counts(&data, 0.001, &energies, 1000.0);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            300.0,
            None,
            vec![0.0005],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_fit_temperature(true)
        .with_transmission_background(BackgroundConfig::default());

        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged);
        let unc = result
            .uncertainties
            .as_ref()
            .expect("KL+bg fit should return density uncertainties");
        assert!(
            unc[0].is_finite() && unc[0] > 0.0,
            "density unc = {}",
            unc[0]
        );
    }

    #[test]
    fn test_lm_uncertainty_not_regressed() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission(&data, 0.001, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.0005],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged);
        let unc = result
            .uncertainties
            .as_ref()
            .expect("LM should still return uncertainties");
        assert!(unc[0].is_finite() && unc[0] > 0.0);
    }

    // ── Energy-scale fitting tests ──

    /// fit_energy_scale + fit_temperature must be rejected.
    #[test]
    fn test_energy_scale_rejects_temperature() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            293.6,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
        .with_fit_temperature(true)
        .with_energy_scale(0.0, 1.0, 25.0);

        let input = InputData::Transmission {
            transmission: t_obs,
            uncertainty: sigma,
        };
        let err = fit_spectrum_typed(&input, &config);
        assert!(err.is_err(), "energy_scale + temperature should fail");
        let msg = err.unwrap_err().to_string();
        assert!(
            msg.contains("fit_energy_scale") && msg.contains("fit_temperature"),
            "error should mention both flags: {msg}"
        );
    }

    /// Energy-scale fitting returns t0_us and l_scale in the result.
    #[test]
    fn test_energy_scale_returns_fitted_params() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            293.6,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
        .with_transmission_background(BackgroundConfig::default())
        .with_energy_scale(0.0, 1.0, 25.0);

        let input = InputData::Transmission {
            transmission: t_obs,
            uncertainty: sigma,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(result.converged, "Fit should converge");
        assert!(
            result.t0_us.is_some(),
            "t0_us should be Some when energy-scale is fitted"
        );
        assert!(
            result.l_scale.is_some(),
            "l_scale should be Some when energy-scale is fitted"
        );
        let t0 = result.t0_us.unwrap();
        let ls = result.l_scale.unwrap();
        // Values should be finite (not NaN/Inf) and within bounds
        assert!(t0.is_finite(), "t0 should be finite, got {t0}");
        assert!(ls.is_finite(), "l_scale should be finite, got {ls}");
        assert!(t0.abs() < 10.0, "t0 should be within bounds, got {t0}");
        assert!(
            ls > 0.98 && ls < 1.02,
            "l_scale should be within bounds, got {ls}"
        );
    }

    /// Without energy-scale fitting, t0_us and l_scale should be None.
    #[test]
    fn test_no_energy_scale_returns_none() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t_obs, sigma) = synthetic_transmission(&data, 0.002, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            293.6,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(Default::default()))
        .with_transmission_background(BackgroundConfig::default());

        let input = InputData::Transmission {
            transmission: t_obs,
            uncertainty: sigma,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert!(
            result.t0_us.is_none(),
            "t0_us should be None without energy-scale"
        );
        assert!(
            result.l_scale.is_none(),
            "l_scale should be None without energy-scale"
        );
    }

    // ==================================================================
    // Joint-Poisson solver integration tests (memo 35 §P1/§P2)
    // ==================================================================

    /// End-to-end: joint-Poisson density recovery at c = 5.98 on synthetic
    /// matched-model counts, via `fit_spectrum_typed`.  Verifies that
    /// `SpectrumFitResult.deviance_per_dof` is populated (P1.2) and that
    /// density is recovered to within 5% on a single-resonance spectrum
    /// under expected (noise-free) counts.
    #[test]
    fn test_joint_poisson_density_recovery_c_5_98() {
        let data = u238_single_resonance();
        let true_density = 0.0005_f64;
        let c = 5.98_f64;
        let lam_ob = 1000.0_f64; // expected open-beam counts per bin (O rate)
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, _) = synthetic_transmission(&data, true_density, &energies);

        // Noise-free expectations under joint-Poisson model:
        //   E[O] = lam_ob, E[S] = c · lam_ob · T
        let open_beam_counts: Vec<f64> = vec![lam_ob; energies.len()];
        let sample_counts: Vec<f64> = t.iter().map(|&ti| c * lam_ob * ti).collect();

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_counts_background(CountsBackgroundConfig {
            c,
            ..Default::default()
        });

        let input = InputData::Counts {
            sample_counts,
            open_beam_counts,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();

        // Deviance-based GOF is populated (P1.2).
        let d_per_dof = result
            .deviance_per_dof
            .expect("joint-Poisson solver must populate deviance_per_dof");
        assert!(
            d_per_dof.is_finite() && d_per_dof >= 0.0,
            "deviance_per_dof = {d_per_dof} is not a valid GOF"
        );
        // On noise-free expected counts, D should be very small (approaches
        // zero in the matched-model limit; allow some slack for numerical
        // error in the forward model).
        assert!(
            d_per_dof < 0.5,
            "noise-free expected-counts fit should give D/dof ≈ 0, got {d_per_dof}"
        );
        // Density recovery.
        let rel_bias = (result.densities[0] - true_density) / true_density;
        assert!(
            rel_bias.abs() < 0.05,
            "density bias {rel_bias} > 5%: fitted={} truth={true_density}",
            result.densities[0]
        );
        // Back-compat: reduced_chi_squared mirrors deviance_per_dof.
        assert!((result.reduced_chi_squared - d_per_dof).abs() < 1e-12);
    }

    /// Counts-KL dispatch rejects `fit_alpha_1` / `fit_alpha_2` — the
    /// profile `λ̂` absorbs the global flux scale (alpha_1 redundant);
    /// alpha_2 / B_det wiring is P3-deferred per memo 35 §P3.
    #[test]
    fn test_joint_poisson_rejects_alpha_fit() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
        let open_beam_counts: Vec<f64> = vec![500.0; energies.len()];
        let sample_counts: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_counts_background(CountsBackgroundConfig {
            fit_alpha_1: true,
            c: 1.0,
            ..Default::default()
        });

        let input = InputData::Counts {
            sample_counts,
            open_beam_counts,
        };
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("fit_alpha_1") || msg.contains("alpha_1"),
            "expected alpha_1 rejection message, got: {msg}"
        );
    }

    /// Transmission + PoissonKL is a valid combination (routes to
    /// `fit_transmission_poisson`, unchanged by the counts-KL collapse).
    /// This test asserts the transmission-KL path still works end-to-end.
    #[test]
    fn test_transmission_poisson_kl_dispatches_to_transmission_path() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, u) = synthetic_transmission(&data, 0.0005, &energies);

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()));

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: u,
        };
        let r = fit_spectrum_typed(&input, &config).unwrap();
        // Transmission-KL path does NOT report deviance_per_dof (that's
        // the counts-domain joint-Poisson GOF only); reduced_chi_squared
        // is the transmission Poisson NLL measure.
        assert!(r.deviance_per_dof.is_none());
        assert!(r.reduced_chi_squared.is_finite());
    }

    // ──────────────────────────────────────────────────────────────────
    // P2.2: transmission_background through the joint-Poisson path.
    // ──────────────────────────────────────────────────────────────────

    /// End-to-end: joint-Poisson with A_n + B_A + B_B + B_C free on
    /// noise-free synthetic counts with known background.  On 201 bins
    /// with 5 free params the (n, A_n) correlation is non-trivial so we
    /// assert the *wiring* is correct (bg reaches the fit, D/dof → 0,
    /// A_n + B_A near truth, density within 10%) rather than EG2-grade
    #[test]
    fn test_joint_poisson_with_transmission_background() {
        let data = u238_single_resonance();
        let true_density = 0.0005_f64;
        let true_anorm = 0.85_f64;
        let true_ba = 0.03_f64;
        let true_bb = -0.01_f64;
        let true_bc = 0.0_f64;
        let c = 5.98_f64;
        let lam_ob = 2000.0_f64;
        let energies: Vec<f64> = (0..201).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t_inner, _) = synthetic_transmission(&data, true_density, &energies);
        let t_out: Vec<f64> = t_inner
            .iter()
            .zip(energies.iter())
            .map(|(&ti, &e)| true_anorm * ti + true_ba + true_bb / e.sqrt() + true_bc * e.sqrt())
            .collect();
        let open_beam_counts: Vec<f64> = vec![lam_ob; energies.len()];
        let sample_counts: Vec<f64> = t_out.iter().map(|&ti| c * lam_ob * ti).collect();

        let bg = BackgroundConfig {
            anorm_init: 1.0,
            back_a_init: 0.0,
            back_b_init: 0.0,
            back_c_init: 0.0,
            back_d_init: 0.01,
            back_f_init: 1.0,
            fit_anorm: true,
            fit_back_a: true,
            fit_back_b: true,
            fit_back_c: true,
            fit_back_d: false,
            fit_back_f: false,
        };

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_counts_background(CountsBackgroundConfig {
            c,
            ..Default::default()
        })
        .with_transmission_background(bg)
        // #486: polish defaults off for real-data regimes where its
        // `fatol = 1e-10` is sub-ULP.  This noise-free synthetic fit
        // is exactly the regime polish was designed for (D → 0
        // achievable, `fatol` physically meaningful) and the test
        // asserts `D/dof < 1.0` which Gauss-Newton alone cannot reach
        // on this 5-free-parameter fit due to the documented n ↔ A_n
        // correlation stall.  Opt in explicitly.
        .with_counts_enable_polish(Some(true));

        let input = InputData::Counts {
            sample_counts,
            open_beam_counts,
        };
        let r = fit_spectrum_typed(&input, &config).unwrap();

        // The invariant P2.2 wiring is supposed to produce is: the 4 bg
        // parameters *actually reach the objective* (the fit produces a
        // near-zero deviance on noise-free expected counts) and the
        // fitter moves them off their initial values.  Density / A_n /
        // B_A recovery at unit-test scale (201 bins, 5 free params)
        // inherits the classic n ↔ A_n correlation — the realistic

        // Deviance-based GOF populated and → 0 on noise-free expected counts.
        // Requires polish (see `with_counts_enable_polish(Some(true))` above).
        let dpd = r.deviance_per_dof.expect("joint-Poisson must report D/dof");
        assert!(
            dpd < 1.0,
            "D/dof = {dpd} unexpectedly large on noise-free fit — bg params not reaching objective?"
        );
        // Density didn't rail to zero.
        assert!(r.densities[0] > 1e-5, "density railed: {}", r.densities[0]);
        // A_n moved off its initial 1.0 toward truth 0.85.
        assert!(
            (r.anorm - 1.0).abs() > 0.05,
            "A_n did not move from init 1.0 (fitted={})",
            r.anorm
        );
        // Background triplet moved off zero at least in one component.
        let bg_moved = r.background.iter().any(|v| v.abs() > 1e-4);
        assert!(
            bg_moved,
            "no bg parameter moved from init 0: {:?}",
            r.background
        );
    }

    /// §P2.2 operational rule: `B_B` or `B_C` free → `B_A` must be free too.
    #[test]
    fn test_joint_poisson_p2_2_requires_back_a_when_back_b_enabled() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
        let ob: Vec<f64> = vec![500.0; energies.len()];
        let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();

        let bg = BackgroundConfig {
            // B_B enabled without B_A → must be rejected.
            fit_anorm: true,
            fit_back_a: false,
            fit_back_b: true,
            fit_back_c: false,
            fit_back_d: false,
            fit_back_f: false,
            ..BackgroundConfig::default()
        };

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_counts_background(CountsBackgroundConfig {
            c: 1.0,
            ..Default::default()
        })
        .with_transmission_background(bg);

        let input = InputData::Counts {
            sample_counts: s,
            open_beam_counts: ob,
        };
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("§P2.2") || msg.contains("B_A"),
            "expected §P2.2 rejection message, got: {msg}"
        );
    }

    /// Joint-Poisson rejects BackD/BackF exponential tail (§P4-deferred).
    #[test]
    fn test_joint_poisson_rejects_back_d_f() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, _) = synthetic_transmission(&data, 0.0005, &energies);
        let ob: Vec<f64> = vec![500.0; energies.len()];
        let s: Vec<f64> = t.iter().map(|&ti| 500.0 * ti).collect();

        let bg = BackgroundConfig {
            fit_back_d: true,
            ..BackgroundConfig::default()
        };

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_counts_background(CountsBackgroundConfig {
            c: 1.0,
            ..Default::default()
        })
        .with_transmission_background(bg);

        let input = InputData::Counts {
            sample_counts: s,
            open_beam_counts: ob,
        };
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            err.to_string().contains("BackD") || err.to_string().contains("§P4"),
            "expected BackD/BackF rejection, got: {err}"
        );
    }

    /// Partial BackD/BackF configurations are rejected with a clear error
    /// on the LM transmission path.  Regression for code-review finding:
    /// pre-fix, `append_background_params` allocated a free index for the
    /// enabled tail parameter but `NormalizedTransmissionModel::new`
    /// (4-term wrapper, fall-back when only one of back_d/back_f was
    /// `Some`) ignored it — the parameter sat at its initial value and
    /// fitter reported it as the "fitted" result.
    #[test]
    fn test_lm_transmission_rejects_partial_back_d_only() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, u) = synthetic_transmission(&data, 0.0005, &energies);

        let bg = BackgroundConfig {
            // BackD enabled, BackF disabled — the partial config the
            // pre-fix code silently accepted.
            fit_back_d: true,
            fit_back_f: false,
            ..BackgroundConfig::default()
        };

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
        .with_transmission_background(bg);

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: u,
        };
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("fit_back_d") && msg.contains("fit_back_f"),
            "expected partial-BackD/F rejection mentioning both flags, got: {msg}"
        );
    }

    /// Symmetric case: BackF enabled without BackD — same rejection.
    #[test]
    fn test_lm_transmission_rejects_partial_back_f_only() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, u) = synthetic_transmission(&data, 0.0005, &energies);

        let bg = BackgroundConfig {
            fit_back_d: false,
            fit_back_f: true,
            ..BackgroundConfig::default()
        };

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
        .with_transmission_background(bg);

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: u,
        };
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("fit_back_d") && msg.contains("fit_back_f"),
            "expected partial-BackD/F rejection mentioning both flags, got: {msg}"
        );
    }

    /// Same rule on the transmission Poisson-KL path.
    #[test]
    fn test_transmission_poisson_kl_rejects_partial_back_d_f() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..51).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let (t, u) = synthetic_transmission(&data, 0.0005, &energies);

        let bg = BackgroundConfig {
            fit_back_d: true,
            fit_back_f: false,
            ..BackgroundConfig::default()
        };

        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_transmission_background(bg);

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: u,
        };
        let err = fit_spectrum_typed(&input, &config).unwrap_err();
        assert!(
            err.to_string().contains("fit_back_d"),
            "expected partial-BackD/F rejection on transmission-PoissonKL path"
        );
    }

    // ------------------------------------------------------------------
    // Fit-energy-range provenance + masked equivalence (#514).
    // ------------------------------------------------------------------

    /// `with_fit_energy_range(...)` round-trips through the accessor.
    #[test]
    fn test_unified_fit_config_fit_energy_range_round_trips() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let cfg = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap();
        assert_eq!(cfg.fit_energy_range(), None);

        let cfg = cfg.with_fit_energy_range(Some((5.0, 50.0))).unwrap();
        assert_eq!(cfg.fit_energy_range(), Some((5.0, 50.0)));

        // Setting back to None clears it.
        let cfg = cfg.with_fit_energy_range(None).unwrap();
        assert_eq!(cfg.fit_energy_range(), None);
    }

    /// `with_fit_energy_range` rejects non-finite or reversed bounds
    /// rather than silently producing an empty active-bin mask
    /// downstream.
    #[test]
    fn test_unified_fit_config_fit_energy_range_rejects_invalid() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..21).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let cfg = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap();

        // Reversed range (lo > hi) is rejected.
        let err = cfg
            .clone()
            .with_fit_energy_range(Some((10.0, 5.0)))
            .unwrap_err();
        assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));

        // Empty range (lo == hi) is rejected — `lo < hi` strictly required.
        let err = cfg
            .clone()
            .with_fit_energy_range(Some((5.0, 5.0)))
            .unwrap_err();
        assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));

        // Non-finite bounds are rejected.
        let err = cfg
            .clone()
            .with_fit_energy_range(Some((f64::NAN, 5.0)))
            .unwrap_err();
        assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));

        let err = cfg
            .clone()
            .with_fit_energy_range(Some((5.0, f64::INFINITY)))
            .unwrap_err();
        assert!(matches!(err, FitConfigError::InvalidFitEnergyRange(_)));
    }

    /// LM fit with `fit_energy_range = Some((min, max))` on the full
    /// grid must yield the same density as the same fit run on the
    /// `[min, max]` slice directly when the residual outside the range
    /// is negligible (SAMMY EMIN/EMAX semantics, paper-acceptance test).
    #[test]
    fn test_fit_energy_range_lm_matches_subset_when_outside_negligible() {
        let data = u238_single_resonance();
        let true_density = 0.002;
        // Wide grid: 0.5..10.5 eV in 0.05 eV steps.  The U-238 single
        // resonance sits well inside [4, 8] eV, so a fit restricted to
        // that range should recover the same density as the full-grid
        // fit (residual outside is negligible).
        let energies: Vec<f64> = (0..201).map(|i| 0.5 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission(&data, true_density, &energies);

        let cfg_full = UnifiedFitConfig::new(
            energies.clone(),
            vec![data.clone()],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));

        let cfg_masked = UnifiedFitConfig::new(
            energies.clone(),
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
        .with_fit_energy_range(Some((4.0, 8.0)))
        .unwrap();

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };
        let r_full = fit_spectrum_typed(&input, &cfg_full).unwrap();
        let r_masked = fit_spectrum_typed(&input, &cfg_masked).unwrap();
        assert!(r_full.converged && r_masked.converged);

        let d_full = r_full.densities[0];
        let d_masked = r_masked.densities[0];
        let rel_err = (d_full - d_masked).abs() / d_full.abs();
        assert!(
            rel_err < 0.01,
            "fit_energy_range LM density {d_masked} should match full-grid {d_full} \
             within 1% (got rel_err = {rel_err})"
        );
    }

    /// Transmission + PoissonKL with `fit_energy_range` set must be
    /// hard-rejected at dispatch.  The legacy `poisson_fit` does not
    /// honour the active-bin mask, so silently fitting on the
    /// (margin-extended) grid would bias the result.  Joint-Poisson
    /// (counts) and LM transmission both honour the mask correctly.
    /// Regression for Round-2 review fix #1 (#514).
    #[test]
    fn test_fit_energy_range_rejected_on_transmission_poisson_path() {
        let data = u238_single_resonance();
        let true_density = 0.002;
        let energies: Vec<f64> = (0..201).map(|i| 0.5 + (i as f64) * 0.05).collect();
        let (t, sigma) = synthetic_transmission(&data, true_density, &energies);

        let cfg = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_fit_energy_range(Some((4.0, 8.0)))
        .unwrap();

        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };
        let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("fit_energy_range")
                && (msg.contains("Poisson") || msg.contains("poisson")),
            "expected rejection message mentioning fit_energy_range and the \
             Poisson-KL path; got: {msg}"
        );
    }

    /// LM dispatcher must reject `fit_energy_range` that selects fewer
    /// than 2 active bins on the configured grid with a clear
    /// "range too narrow" error — instead of a confusing non-converged
    /// LM result.  Regression for #517 R3 P2 (Copilot finding #2).
    #[test]
    fn test_fit_energy_range_lm_rejects_too_narrow() {
        let data = u238_single_resonance();
        // 0.5..10.5 eV in 0.5 eV steps.  Range [4.6, 4.7] selects
        // exactly zero bins (no grid point inside; closest are 4.5
        // and 5.0).
        let energies: Vec<f64> = (0..21).map(|i| 0.5 + (i as f64) * 0.5).collect();
        let (t, sigma) = synthetic_transmission(&data, 0.002, &energies);
        let cfg = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
        .with_fit_energy_range(Some((4.6, 4.7)))
        .unwrap();
        let input = InputData::Transmission {
            transmission: t,
            uncertainty: sigma,
        };
        let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("fit_energy_range") && msg.contains("active bin"),
            "expected too-narrow-range rejection; got: {msg}"
        );
    }

    /// Joint-Poisson dispatcher must reject too-narrow ranges with the
    /// same clear error.  Regression for #517 R3 P2 (Copilot #3).
    #[test]
    fn test_fit_energy_range_jp_rejects_too_narrow() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..21).map(|i| 0.5 + (i as f64) * 0.5).collect();
        let (sample, open_beam) = synthetic_counts(&data, 0.002, &energies, 1000.0);
        let cfg = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_fit_energy_range(Some((4.6, 4.7)))
        .unwrap();
        let input = InputData::Counts {
            sample_counts: sample,
            open_beam_counts: open_beam,
        };
        let err = fit_spectrum_typed(&input, &cfg).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("fit_energy_range") && msg.contains("active bin"),
            "expected too-narrow-range rejection; got: {msg}"
        );
    }

    // ── Issue #608 (PR #609): Rust coverage for paths previously exercised only
    //    via Python / the spatial-identity path (codecov/patch gap) ───────────

    /// `evaluate_jacobian_and_fisher` (the research Fisher / `compute_model_jacobian`
    /// entry point) was exercised only through the Python bindings, which the
    /// Rust-only coverage run excludes — so its body, and the #608
    /// compute-working-grid-σ branch, were uncovered.  Drive it with a
    /// Gaussian-resolution config carrying NO precomputed σ: that exercises the
    /// aux-grid σ build + the model construction + the analytical Jacobian /
    /// Fisher assembly.
    #[test]
    fn evaluate_jacobian_and_fisher_gaussian_aux_grid() {
        use nereids_physics::resolution::{ResolutionFunction, ResolutionParams};
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let n_e = energies.len();
        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            300.0,
            Some(ResolutionFunction::Gaussian(
                ResolutionParams::new(25.0, 0.5, 0.005, 0.0).unwrap(),
            )),
            vec![0.001],
        )
        .unwrap();
        let flux = vec![5000.0; n_e];
        let background = vec![10.0; n_e];
        let result = evaluate_jacobian_and_fisher(&config, &flux, &background).unwrap();
        assert_eq!(result.model_prediction.len(), n_e);
        assert!(result.model_prediction.iter().all(|v| v.is_finite()));
        assert_eq!(result.param_names.len(), 1, "one free density parameter");
        // Density Fisher information must be finite + positive (well-posed
        // measurement); the Jacobian endpoints must be finite.
        let f00 = result.fisher.get(0, 0);
        assert!(f00.is_finite() && f00 > 0.0, "Fisher[0,0] = {f00}");
        assert!(result.jacobian.get(0, 0).is_finite());
        assert!(result.jacobian.get(n_e - 1, 0).is_finite());
    }

    /// Non-spatial energy-scale fit through `fit_spectrum_typed` (LM
    /// transmission): exercises `seed_energy_scale_in_params` +
    /// `build_energy_scale_transmission_model` + the energy-scale LM path +
    /// `t0_us` / `l_scale` result population.  Prior Rust coverage was
    /// spatial-only on IDENTITY data; here we inject a NON-identity (t0,
    /// L_scale) and confirm the wiring recovers the density (the (t0, L_scale)
    /// valley is shallow, so we assert the physically-observable density, not
    /// tight individual parameters — cf. the Python `TestFitEnergyScaleRecovery`
    /// rationale).
    #[test]
    fn fit_spectrum_typed_energy_scale_lm_recovers_calibration() {
        let data = hf178_mlbw_two_resonances(); // s-waves @ 7.8 + 16.9 eV
        let energies: Vec<f64> = (0..700).map(|i| 4.0 + (i as f64) * 0.025).collect();
        let true_density = 0.05_f64;
        let (t0_true, ls_true) = (1.5_f64, 1.004_f64);
        // Synthesize measured data with the injected calibration via the model
        // (no resolution: dip POSITIONS drive the seed).
        let model = EnergyScaleTransmissionModel::new(
            Arc::new(vec![data.clone()]),
            Arc::new(vec![0]),
            Arc::new(vec![1.0]),
            293.6,
            energies.clone(),
            25.0,
            1,
            2,
            None,
        );
        let measured = model.evaluate(&[true_density, t0_true, ls_true]).unwrap();
        let uncertainty = vec![1e-3; energies.len()];
        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["Hf-178".into()],
            293.6,
            None,
            vec![0.04], // start density away from truth (0.05) so recovery is real
        )
        .unwrap()
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()))
        .with_energy_scale(0.0, 1.0, 25.0); // cold (t0,L) start ⇒ the seed must correct it
        let input = InputData::Transmission {
            transmission: measured,
            uncertainty,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        let t0 = result
            .t0_us
            .expect("t0_us populated when fit_energy_scale=true");
        let ls = result
            .l_scale
            .expect("l_scale populated when fit_energy_scale=true");
        assert!(t0.is_finite() && ls.is_finite(), "t0={t0}, L={ls}");
        assert!(result.converged, "energy-scale LM fit should converge");
        assert!(
            (result.densities[0] - true_density).abs() / true_density < 0.10,
            "density: fitted={}, true={true_density}",
            result.densities[0]
        );
    }

    /// Non-spatial counts-KL energy-scale fit through `fit_spectrum_typed`:
    /// exercises `seed_energy_scale_in_params`' KL transmission proxy
    /// `(sample − bg)/(flux − bg)` + the counts-KL energy-scale dispatch —
    /// uncovered by the Rust suite (Python + spatial-identity only).
    #[test]
    fn fit_spectrum_typed_energy_scale_counts_kl_seeds_via_proxy() {
        let data = hf178_mlbw_two_resonances();
        let energies: Vec<f64> = (0..700).map(|i| 4.0 + (i as f64) * 0.025).collect();
        let true_density = 0.05_f64;
        let (t0_true, ls_true) = (1.0_f64, 1.003_f64);
        let model = EnergyScaleTransmissionModel::new(
            Arc::new(vec![data.clone()]),
            Arc::new(vec![0]),
            Arc::new(vec![1.0]),
            293.6,
            energies.clone(),
            25.0,
            1,
            2,
            None,
        );
        let t = model.evaluate(&[true_density, t0_true, ls_true]).unwrap();
        // Counts: sample = flux·T.  Zero detector background — the joint-Poisson
        // KL path does not yet wire B_det (memo 35 §P3.2); the seed proxy
        // (sample − 0)/(flux − 0) = T still reconstructs the dip positions.
        let flux: Vec<f64> = vec![5000.0; energies.len()];
        let background: Vec<f64> = vec![0.0; energies.len()];
        let sample: Vec<f64> = t
            .iter()
            .zip(flux.iter())
            .map(|(&ti, &fi)| fi * ti)
            .collect();
        let config = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["Hf-178".into()],
            293.6,
            None,
            vec![0.04],
        )
        .unwrap()
        .with_solver(SolverConfig::PoissonKL(PoissonConfig::default()))
        .with_energy_scale(0.0, 1.0, 25.0);
        let input = InputData::CountsWithNuisance {
            sample_counts: sample,
            flux,
            background,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        let t0 = result
            .t0_us
            .expect("t0_us populated when fit_energy_scale=true");
        let ls = result
            .l_scale
            .expect("l_scale populated when fit_energy_scale=true");
        assert!(t0.is_finite() && ls.is_finite(), "t0={t0}, L={ls}");
    }

    /// The #608 working-grid σ + layout validation in
    /// `validate_precomputed_cross_sections` (the up-front guard for the
    /// Gaussian aux-grid path) — every malformed-input error branch, so a bad
    /// `with_precomputed_work_cross_sections` setter call surfaces a typed
    /// `ShapeMismatch` instead of a per-pixel panic.
    #[test]
    fn validate_precomputed_work_cross_sections_error_branches() {
        use nereids_physics::transmission::WorkingGridLayout;
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..11).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let n_e = energies.len();
        let n_work = n_e + 2; // a non-identity aux grid
        let good_layout = || {
            Arc::new(WorkingGridLayout {
                energies: (0..n_work).map(|i| 1.0 + (i as f64) * 0.09).collect(),
                data_indices: (0..n_e).collect(),
            })
        };
        let cfg = |work_xs: Vec<Vec<f64>>, layout: Arc<WorkingGridLayout>| {
            UnifiedFitConfig::new(
                energies.clone(),
                vec![data.clone()],
                vec!["U-238".into()],
                0.0,
                None,
                vec![0.001],
            )
            .unwrap()
            .with_precomputed_cross_sections(Arc::new(vec![vec![1.0f64; n_e]])) // valid data-grid σ
            .with_precomputed_work_cross_sections(Arc::new(work_xs), layout)
        };
        let expect =
            |c: &UnifiedFitConfig, needle: &str| match validate_precomputed_cross_sections(c) {
                Err(PipelineError::ShapeMismatch(m)) => {
                    assert!(m.contains(needle), "expected {needle:?}, got: {m}")
                }
                other => panic!("expected ShapeMismatch({needle:?}), got {other:?}"),
            };
        // empty working σ
        expect(&cfg(vec![], good_layout()), "must not be empty");
        // working-σ row length != working-grid length
        expect(
            &cfg(vec![vec![1.0; n_work - 1]], good_layout()),
            "working-grid energies",
        );
        // non-finite working σ
        let mut nan_row = vec![1.0; n_work];
        nan_row[3] = f64::NAN;
        expect(&cfg(vec![nan_row], good_layout()), "non-finite");
        // working-σ row count != data-grid σ row count (2 work rows vs 1 data row)
        expect(
            &cfg(vec![vec![1.0; n_work], vec![1.0; n_work]], good_layout()),
            "rows but",
        );
        // layout maps a different number of data points than config.energies
        let short = Arc::new(WorkingGridLayout {
            energies: (0..n_work).map(|i| 1.0 + (i as f64) * 0.09).collect(),
            data_indices: (0..n_e - 1).collect(),
        });
        expect(&cfg(vec![vec![1.0; n_work]], short), "layout maps");
        // layout index out of range for the working grid
        let oor = Arc::new(WorkingGridLayout {
            energies: (0..n_work).map(|i| 1.0 + (i as f64) * 0.09).collect(),
            data_indices: {
                let mut v: Vec<usize> = (0..n_e).collect();
                v[0] = n_work + 5;
                v
            },
        });
        expect(&cfg(vec![vec![1.0; n_work]], oor), "out of");
    }

    /// Cover the OTHER branches of the #608 `evaluate_jacobian_and_fisher` σ
    /// restructure: the identity-layout path (no resolution ⇒ working grid ==
    /// data grid) and the early-return when `fit_temperature` is set (the
    /// `TransmissionFitModel` builds its own working-grid base σ).
    #[test]
    fn evaluate_jacobian_and_fisher_identity_and_temperature_paths() {
        let data = u238_single_resonance();
        let energies: Vec<f64> = (0..101).map(|i| 1.0 + (i as f64) * 0.1).collect();
        let n_e = energies.len();
        let flux = vec![5000.0; n_e];
        let background = vec![10.0; n_e];
        // (a) No resolution ⇒ identity working-grid layout.
        let cfg_id = UnifiedFitConfig::new(
            energies.clone(),
            vec![data.clone()],
            vec!["U-238".into()],
            300.0,
            None,
            vec![0.001],
        )
        .unwrap();
        let r_id = evaluate_jacobian_and_fisher(&cfg_id, &flux, &background).unwrap();
        assert_eq!(r_id.model_prediction.len(), n_e);
        let f00 = r_id.fisher.get(0, 0);
        assert!(
            f00.is_finite() && f00 > 0.0,
            "identity-path Fisher[0,0]={f00}"
        );
        // (b) fit_temperature ⇒ the σ-branch early-returns `config`;
        //     TransmissionFitModel builds its own working-grid base σ.
        let cfg_t = UnifiedFitConfig::new(
            energies,
            vec![data],
            vec!["U-238".into()],
            300.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_fit_temperature(true);
        let r_t = evaluate_jacobian_and_fisher(&cfg_t, &flux, &background).unwrap();
        assert_eq!(
            r_t.param_names.len(),
            2,
            "density + temperature free params"
        );
        assert!(r_t.model_prediction.iter().all(|v| v.is_finite()));
    }

    /// Cover `build_transmission_model`'s group-collapse path: a grouped config
    /// (two isotopes → one density param) carrying PRECOMPUTED per-member σ is
    /// collapsed to per-group effective σ before the model build (the grouped +
    /// precomputed path the spatial pipeline drives per pixel).
    #[test]
    fn fit_spectrum_typed_grouped_precomputed_collapses_by_groups() {
        use nereids_core::types::{Isotope, IsotopeGroup};
        let rd1 = synthetic_single_resonance(92, 235, 233.025, 5.0);
        let rd2 = synthetic_single_resonance(92, 238, 236.006, 7.0);
        let energies: Vec<f64> = (0..301).map(|i| 1.0 + (i as f64) * 0.05).collect();
        let true_density = 0.0005_f64;
        let sample = phys_transmission::SampleParams::new(
            0.0,
            vec![
                (rd1.clone(), true_density * 0.6),
                (rd2.clone(), true_density * 0.4),
            ],
        )
        .unwrap();
        let transmission = phys_transmission::forward_model(&energies, &sample, None).unwrap();
        let uncertainty: Vec<f64> = transmission.iter().map(|&t| 0.01 * t.max(0.01)).collect();
        // Per-member precomputed σ (2 rows) ⇒ triggers the group collapse.
        let per_member = phys_transmission::broadened_cross_sections(
            &energies,
            &[rd1.clone(), rd2.clone()],
            0.0,
            None,
            None,
        )
        .unwrap();
        let iso1 = Isotope::new(92, 235).unwrap();
        let iso2 = Isotope::new(92, 238).unwrap();
        let group = IsotopeGroup::custom("U".into(), vec![(iso1, 0.6), (iso2, 0.4)]).unwrap();
        let config = UnifiedFitConfig::new(
            energies,
            vec![rd1.clone()],
            vec!["placeholder".into()],
            0.0,
            None,
            vec![0.001],
        )
        .unwrap()
        .with_groups(&[(&group, &[rd1, rd2])], vec![0.001])
        .unwrap()
        .with_precomputed_cross_sections(Arc::new(per_member))
        .with_solver(SolverConfig::LevenbergMarquardt(LmConfig::default()));
        let input = InputData::Transmission {
            transmission,
            uncertainty,
        };
        let result = fit_spectrum_typed(&input, &config).unwrap();
        assert_eq!(result.densities.len(), 1, "one group density");
        assert!(
            (result.densities[0] - true_density).abs() / true_density < 0.01,
            "grouped+precomputed density: fitted={}, true={true_density}",
            result.densities[0]
        );
    }

    /// `peak_match_energy_scale_seed` rejects (→ keeps the cold start) a fit that
    /// lands OUTSIDE the parameter bounds rather than clamping onto a bound
    /// (issue #608 R2): inject an L_scale (1.03) beyond the seed's (0.99, 1.01).
    #[test]
    fn peak_match_energy_scale_seed_rejects_out_of_bounds() {
        let data = hf178_mlbw_two_resonances();
        let energies: Vec<f64> = (0..900).map(|i| 4.0 + (i as f64) * 0.02).collect();
        let model = EnergyScaleTransmissionModel::new(
            Arc::new(vec![data.clone()]),
            Arc::new(vec![0]),
            Arc::new(vec![1.0]),
            293.6,
            energies.clone(),
            25.0,
            1,
            2,
            None,
        );
        // L_scale = 1.03 is well outside (0.99, 1.01); the seed would fit ≈1.03.
        let t_obs = model.evaluate(&[0.05, 0.0, 1.03]).unwrap();
        let config = UnifiedFitConfig::new(
            energies.clone(),
            vec![data],
            vec!["Hf-178".into()],
            293.6,
            None,
            vec![0.05],
        )
        .unwrap()
        .with_energy_scale(0.0, 1.0, 25.0);
        assert!(
            peak_match_energy_scale_seed(
                &t_obs,
                config.energies(),
                &config,
                25.0,
                (-10.0, 10.0),
                (0.99, 1.01),
            )
            .is_none(),
            "an out-of-bounds calibration fit must be rejected (cold-start fallback)"
        );
    }
}