rexafs 0.2.5

#![allow(dead_code)]
#![allow(unused_imports)]
#![allow(unused_variables)]

// Import standard library dependencies
#[cfg(feature = "trust-region")]
use std::collections::HashMap;
use std::collections::{hash_map::DefaultHasher, VecDeque};
use std::error::Error;
use std::hash::Hasher;
use std::ops::Deref;
use std::sync::{Arc, Mutex, OnceLock};

// Import external dependencies
use super::spline;
#[cfg(feature = "trust-region")]
use apex_solver::core::problem::Problem as ApexProblem;
#[cfg(feature = "trust-region")]
use apex_solver::factors::Factor as ApexFactor;
#[cfg(feature = "trust-region")]
use apex_solver::linalg::{JacobianMode as ApexJacobianMode, LinearSolverType};
#[cfg(feature = "trust-region")]
use apex_solver::manifold::ManifoldType;
#[cfg(feature = "trust-region")]
use apex_solver::optimizer::dog_leg::{DogLeg, DogLegConfig};
use levenberg_marquardt::{LeastSquaresProblem, LevenbergMarquardt};
use nalgebra::{DMatrix, DVector, Dyn, Owned};
#[cfg(feature = "trust-region")]
use nalgebra_apex::DVector as ApexDVector;
use serde::{Deserialize, Serialize};

// Import internal dependencies
use super::errors::{BackgroundError, DataError};
use super::lmutils::LMParameters;
use super::mathutils::{self, splev_jacobian, MathUtils};
use super::normalization::{self, Normalization};
use super::xafsutils::FTWindow;
use super::xrayfft::{FFTUtils, XFFTReverse, XFFT};
use super::{xafsutils, xrayfft};

mod fixed;

/// Weak, constant endpoint penalty for the mean-square fixed AUTOBK objective.
pub const DEFAULT_CLAMP_LAMBDA: f64 = 0.001;

const DEFAULT_LINEAR_REGULARIZATION: f64 = 1.0e-4;
const DEFAULT_LINEAR_CONDITION_LIMIT: f64 = 1.0e8;
const DEFAULT_LINEAR_RESIDUAL_RATIO_LIMIT: f64 = 1.05;
const AUTOBK_WORKSPACE_CACHE_CAPACITY: usize = 8;

/// Solver backend for the spline coefficients.
///
/// LinearDirect with FixedPenalty is the recommended default; iterative solvers
/// require a legacy clamp policy and optimize a different residual.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub enum AUTOBKSolver {
    /// Legacy iterative trust-region dog-leg; requires the `trust-region` feature.
    TrustRegionDogLeg,
    /// Legacy Levenberg–Marquardt iteration with the dynamic clamp residual.
    LegacyLm,
    /// Solve a linear system; recommended with the fixed endpoint objective.
    LinearDirect,
}

#[cfg(feature = "trust-region")]
fn default_autobk_fallback_solver() -> AUTOBKSolver {
    AUTOBKSolver::TrustRegionDogLeg
}

#[cfg(not(feature = "trust-region"))]
fn default_autobk_fallback_solver() -> AUTOBKSolver {
    AUTOBKSolver::LegacyLm
}

/// Endpoint clamp convention. The fixed-penalty default and historical dynamic
/// models have different objectives and should not be interpreted as solver-only choices.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
pub enum AUTOBKClampScalePolicy {
    /// One linear solve of mean(low-R residual²) + lambda * mean(active endpoint χ²).
    #[default]
    FixedPenalty,
    /// Legacy scale frozen from the initial spline, with ridge retries/fallback.
    Fixed,
    /// Legacy direct solve, refresh the residual-dependent scale, then solve again.
    TwoPass,
}

fn legacy_clamp_policy() -> Option<AUTOBKClampScalePolicy> {
    Some(AUTOBKClampScalePolicy::Fixed)
}

#[derive(Debug, Clone)]
struct AutobkLinearWorkspace {
    signature: u64,
    design_matrix: Arc<DMatrix<f64>>,
}

fn autobk_workspace_cache() -> &'static Mutex<VecDeque<AutobkLinearWorkspace>> {
    static CACHE: OnceLock<Mutex<VecDeque<AutobkLinearWorkspace>>> = OnceLock::new();
    CACHE.get_or_init(|| Mutex::new(VecDeque::new()))
}

/// Background algorithm and its settings/results.
///
/// AUTOBK implements spline-based low-R suppression as introduced by
/// [Newville et al. (1993)](https://doi.org/10.1103/PhysRevB.47.14126).
/// ILPBkg is an unimplemented placeholder; None explicitly skips calculation.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
#[allow(clippy::large_enum_variant)]
pub enum BackgroundMethod {
    /// Implemented spline-based background removal.
    AUTOBK(AUTOBK),
    /// Retained placeholder; background calculation returns `NotImplemented`.
    ILPBkg(ILPBkg),
    /// Explicitly skip background calculation; no k or chi arrays are produced.
    None,
}

impl From<AUTOBK> for BackgroundMethod {
    fn from(parameters: AUTOBK) -> Self {
        Self::AUTOBK(parameters)
    }
}

// Allow direct settings while preserving existing optional-enum setter calls.
impl From<AUTOBK> for Option<BackgroundMethod> {
    fn from(parameters: AUTOBK) -> Self {
        Some(parameters.into())
    }
}

impl Default for BackgroundMethod {
    fn default() -> Self {
        BackgroundMethod::AUTOBK(AUTOBK::default())
    }
}

impl BackgroundMethod {
    /// Select AUTOBK with recommended settings.
    pub fn new() -> BackgroundMethod {
        BackgroundMethod::AUTOBK(AUTOBK::new())
    }

    /// Select AUTOBK with recommended settings.
    pub fn new_autobk() -> BackgroundMethod {
        BackgroundMethod::AUTOBK(AUTOBK::new())
    }

    /// Create the ILPBkg placeholder; background calculation returns `NotImplemented`.
    pub fn new_ilpbkg() -> BackgroundMethod {
        BackgroundMethod::ILPBkg(ILPBkg::new())
    }

    /// Fit a smooth background from borrowed energy (eV) and absorption arrays.
    ///
    /// Calculates missing normalization, mutates the supplied normalization settings
    /// and replaces this object's cached outputs. Inputs must be finite, matched
    /// and increasing with enough post-edge points. Invalid cutoffs, spline geometry
    /// or failed solves return `BackgroundError`; unavailable algorithms return
    /// `NotImplemented`. Input arrays themselves are unchanged.
    pub fn calc_background(
        &mut self,
        energy: &DVector<f64>,
        mu: &DVector<f64>,
        normalization_param: &mut Option<normalization::NormalizationMethod>,
    ) -> Result<&mut Self, BackgroundError> {
        match self {
            BackgroundMethod::AUTOBK(autobk) => {
                autobk.calc_background(energy, mu, normalization_param)?;
                Ok(self)
            }
            BackgroundMethod::ILPBkg(_ilpbkg) => Err(BackgroundError::NotImplemented {
                feature: "ILPBkg background removal".to_string(),
            }),
            BackgroundMethod::None => Ok(self),
        }
    }

    /// Return a newly allocated copy without recalculating.
    ///
    /// Uniform output wave-number grid in Å⁻¹ starting at zero.
    /// `None` before calculation.
    pub fn get_k(&self) -> Option<DVector<f64>> {
        match self {
            BackgroundMethod::AUTOBK(autobk) => autobk.get_k().cloned(),
            _ => None,
        }
    }

    /// Return a newly allocated copy without recalculating.
    ///
    /// Dimensionless unweighted EXAFS on `k`, without objective-only standard
    /// subtraction. `None` before calculation.
    pub fn get_chi(&self) -> Option<DVector<f64>> {
        match self {
            BackgroundMethod::AUTOBK(autobk) => autobk.get_chi().cloned(),
            _ => None,
        }
    }
}

/// Cubic-spline absorption background removal by suppressing low-R EXAFS.
///
/// Start with `new()` and the recommended `rbkg=1.0` Å, `kstep=0.05` Å⁻¹,
/// LinearDirect solver and FixedPenalty endpoint strength 0.001. The original
/// AUTOBK method is described by
/// [Newville et al. (1993)](https://doi.org/10.1103/PhysRevB.47.14126).
/// FixedPenalty is a rexafs modification: one linear least-squares solve of
/// mean squared low-R residual plus a fixed weighted endpoint penalty.
///
/// For fixed grid/window/edge step, chi(c)=y−B*c is affine in cubic spline
/// coefficients c. Here y is normalized data, optionally minus an objective-only
/// standard, and B contains the spline basis on the output k grid. If h(c)
/// stacks m real/imaginary low-R Fourier entries, the objective is
///
/// ```text
/// J(c) = sum(h(c)^2)/m
///      + lambda * sum(active endpoint weight^2 * chi(c)^2)/N_active
/// ```
///
/// N_active counts enabled endpoint rows; zero weight omits an endpoint and
/// lambda=0 or nclamp=0 omits the penalty. Each enabled end contributes up to
/// nclamp output samples. A row may appear at both ends when the sets overlap;
/// those occurrences are counted separately. lambda is a numerical penalty
/// strength tied to the chosen Fourier/k-weight units, not a universal
/// dimensionless physical quantity or confidence level. The fixed residual
/// FFT uses 0.05/sqrt(pi), while legacy policies use the actual kstep/sqrt(pi).
///
/// The single least-squares system stacks h and endpoint rows multiplied by
/// sqrt(lambda*m/N_active). Column scaling and singular-value decomposition
/// solve it without adding regularization. Rank, conditioning and stationarity
/// failures return errors instead of changing lambda or falling back.
/// This fixed objective is implemented in the private `background::fixed` module.
///
/// `calc_background` computes missing normalization and stores new background,
/// k and chi arrays. Directly editing fields does not invalidate existing
/// spectrum results; use `Spectrum::set_background_method` for managed updates.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
#[serde(default)]
pub struct AUTOBK {
    /// Background energy origin in eV. `None` uses the normalization E0. An explicit
    /// value changes conversion to k without replacing a valid normalization E0.
    pub ek0: Option<f64>,
    /// Low-R background cutoff in Å; recommended default 1.0, finite and positive.
    /// Larger values increase spline flexibility and may remove structural signal.
    /// This is an algorithmic cutoff, not an estimate of a bond length.
    pub rbkg: Option<f64>,
    /// Requested number of fitted cubic-spline coefficients. `None` uses
    /// `1 + floor(2*rbkg*(kmax-kmin)/pi)`; the final count is bounded to 5–128.
    /// Changing it overrides flexibility but does not independently change the low-R cutoff.
    pub nknots: Option<i32>,
    /// Lower spline/window boundary in Å⁻¹; default 0.0. The returned k grid still
    /// starts at zero; endpoint penalties use that output grid, not this boundary.
    pub kmin: Option<f64>,
    /// Upper spline/window boundary in Å⁻¹. `None` uses the available post-edge range;
    /// an explicit value is capped at that range. Must exceed `kmin`.
    pub kmax: Option<f64>,
    /// Uniform output k spacing in Å⁻¹; default 0.05, finite and positive.
    /// This sets the physical R grid. FixedPenalty uses the fixed numerical
    /// amplitude reference `0.05 / sqrt(pi)`; legacy policies instead use the
    /// actual `kstep / sqrt(pi)`.
    pub kstep: Option<f64>,
    /// Number of output k samples penalized at each enabled end; default 3.
    /// FixedPenalty uses up to the available point count and includes the last sample.
    /// Zero disables endpoint rows; FixedPenalty rejects negative counts.
    pub nclamp: Option<i32>,
    /// Low-k endpoint weight; default 0 (disabled). FixedPenalty uses its absolute
    /// value and squares that weight in the objective, so doubling it quadruples
    /// that end's contribution.
    pub clamp_lo: Option<i32>,
    /// High-k endpoint weight; default 1. FixedPenalty uses its absolute value;
    /// zero disables that end. The returned chi is not multiplied by this weight.
    pub clamp_hi: Option<i32>,
    /// Internal Fourier length; default 2048. Use enough points for the complete
    /// output grid; the low-level FFT pads short arrays and truncates longer arrays.
    pub nfft: Option<i32>,
    /// Objective-only standard, interpolated from `k_std` onto the output k grid.
    /// FixedPenalty expects dimensionless chi. Legacy policies subtract the
    /// standard before edge-step division, so they require values in original
    /// mu units. It is never subtracted from the returned chi.
    pub chi_std: Option<DVector<f64>>,
    /// Wave-number grid of `chi_std`, in Å⁻¹; both arrays must be supplied together.
    /// Use finite matching arrays with strictly increasing k for interpolation.
    pub k_std: Option<DVector<f64>>,
    /// Integer k-weight used in the background objective; default 1.
    /// It is separate from public forward-transform weighting and does not multiply
    /// the returned unweighted chi. Larger weights emphasize high-k data and noise.
    pub kweight: Option<i32>,
    /// Window family for the background objective; default Hanning.
    /// Changing the family alters the objective and can affect the fitted background.
    pub window: FTWindow,
    /// Internal window parameter; default 0.1, applied at both ends.
    /// Usually a taper width in Å⁻¹; Kaiser–Bessel also uses its value for shape,
    /// Gaussian uses it for width, and fractional Hanning uses fractional geometry.
    pub dk: Option<f64>,
    /// Spline solver; recommended default `LinearDirect`. FixedPenalty requires
    /// this solver and returns failures as errors instead of falling back.
    pub solver: Option<AUTOBKSolver>,
    #[serde(default = "legacy_clamp_policy")]
    /// Endpoint objective; recommended default `FixedPenalty`. `Fixed` and
    /// `TwoPass` retain legacy scale models. Older serialized settings lacking this
    /// field restore `Fixed` for historical compatibility.
    pub clamp_scale_policy: Option<AUTOBKClampScalePolicy>,
    /// FixedPenalty endpoint strength; recommended default 0.001. Must be finite
    /// and nonnegative; zero removes endpoint rows. Its numerical balance is tied
    /// to the fixed 0.05/sqrt(pi) FFT scale, k-weight and window. This is an
    /// empirical rexafs choice, not a universal dimensionless physical constant
    /// or the dynamic clamp model in stock Larch.
    pub clamp_lambda: Option<f64>,
    /// Legacy Fixed/TwoPass ridge magnitude; default 1e-4.
    /// Unused by FixedPenalty, which never adds a ridge to force a solution.
    pub linear_regularization: Option<f64>,
    /// Direct-solver condition threshold; default 1e8. FixedPenalty checks the
    /// ratio of largest to smallest singular value after column scaling and requires
    /// a finite limit at least 1, including when reusing a cached factorization.
    pub linear_condition_limit: Option<f64>,
    /// Legacy Fixed/TwoPass acceptance limit for solved/base residual norm;
    /// default 1.05. Unused by FixedPenalty.
    pub linear_residual_ratio_limit: Option<f64>,
    /// Allow legacy direct-solve failures to use `linear_fallback_solver`; default
    /// true. Despite the historical name, the chosen fallback may be dog-leg.
    /// Unused by FixedPenalty, which returns failed solves as errors.
    pub linear_fallback_to_lm: Option<bool>,
    /// Legacy nonlinear fallback. New defaults use TrustRegionDogLeg with the
    /// `trust-region` feature and LegacyLm otherwise. When this field is explicitly
    /// `None`, legacy fill behavior selects LegacyLm unless fallback is disabled.
    /// Unused by FixedPenalty.
    pub linear_fallback_solver: Option<AUTOBKSolver>,
    /// Reuse compatible direct-solver geometry; default true. FixedPenalty also
    /// caches column scaling and SVD factors, but always solves a new right-hand side
    /// for each spectrum. This does not reuse another spectrum's fitted answer.
    pub linear_workspace_cache: Option<bool>,
    /// Smooth background on the original energy grid, in input mu units.
    /// Outside the fitted interval this equals the original mu; `None` before calculation.
    pub bkg: Option<DVector<f64>>,
    /// Dimensionless `(mu - bkg) / edge_step` on the original energy grid.
    /// It is zero outside the fitted interval; `None` before calculation.
    pub chie: Option<DVector<f64>>,
    /// Uniform output wave-number grid in Å⁻¹ starting at zero.
    /// `None` before calculation.
    pub k: Option<DVector<f64>>,
    /// Dimensionless unweighted EXAFS on `k`, without objective-only standard
    /// subtraction. `None` before calculation.
    pub chi: Option<DVector<f64>>,
}

impl Default for AUTOBK {
    fn default() -> Self {
        AUTOBK {
            ek0: None,
            rbkg: Some(1.0),
            nknots: None,
            kmin: Some(0.0),
            kmax: None,
            kstep: Some(0.05),
            nclamp: Some(3),
            clamp_lo: Some(0),
            clamp_hi: Some(1),
            nfft: Some(2048),
            chi_std: None,
            k_std: None,
            kweight: Some(1),
            window: FTWindow::Hanning,
            dk: Some(0.1),
            solver: Some(AUTOBKSolver::LinearDirect),
            clamp_scale_policy: Some(AUTOBKClampScalePolicy::FixedPenalty),
            clamp_lambda: Some(DEFAULT_CLAMP_LAMBDA),
            linear_regularization: Some(DEFAULT_LINEAR_REGULARIZATION),
            linear_condition_limit: Some(DEFAULT_LINEAR_CONDITION_LIMIT),
            linear_residual_ratio_limit: Some(DEFAULT_LINEAR_RESIDUAL_RATIO_LIMIT),
            linear_fallback_to_lm: Some(true),
            linear_fallback_solver: Some(default_autobk_fallback_solver()),
            linear_workspace_cache: Some(true),
            bkg: None,
            chie: None,
            k: None,
            chi: None,
        }
    }
}

impl AUTOBK {
    /// Create recommended background settings; equivalent to `default()`.
    pub fn new() -> AUTOBK {
        AUTOBK::default()
    }

    /// Resolve unset scalar defaults without fitting a background.
    /// Data-dependent E0, maximum k and coefficient count are resolved later.
    /// Filled scalar values are retained. calc_background also stores its
    /// resolved ek0, but automatic kmax/nknots remain None and are calculated
    /// locally from the current data.
    pub fn fill_parameter(&mut self) -> Result<(), BackgroundError> {
        if self.rbkg.is_none() {
            self.rbkg = Some(1.0);
        }

        if self.kmin.is_none() {
            self.kmin = Some(0.0);
        }

        if self.kstep.is_none() {
            self.kstep = Some(0.05);
        }

        if self.nclamp.is_none() {
            self.nclamp = Some(3);
        }

        if self.clamp_lo.is_none() {
            self.clamp_lo = Some(0);
        }

        if self.clamp_hi.is_none() {
            self.clamp_hi = Some(1);
        }

        if self.nfft.is_none() {
            self.nfft = Some(2048);
        }

        if self.kweight.is_none() {
            self.kweight = Some(1);
        }

        if self.dk.is_none() {
            self.dk = Some(0.1);
        }

        if self.solver.is_none() {
            self.solver = Some(AUTOBKSolver::LinearDirect);
        }

        if self.clamp_scale_policy.is_none() {
            self.clamp_scale_policy = Some(AUTOBKClampScalePolicy::FixedPenalty);
        }

        if self.clamp_lambda.is_none() {
            self.clamp_lambda = Some(DEFAULT_CLAMP_LAMBDA);
        }

        if self.linear_regularization.is_none() {
            self.linear_regularization = Some(DEFAULT_LINEAR_REGULARIZATION);
        }

        if self.linear_condition_limit.is_none() {
            self.linear_condition_limit = Some(DEFAULT_LINEAR_CONDITION_LIMIT);
        }

        if self.linear_residual_ratio_limit.is_none() {
            self.linear_residual_ratio_limit = Some(DEFAULT_LINEAR_RESIDUAL_RATIO_LIMIT);
        }

        if self.linear_fallback_solver.is_none() {
            let fallback_solver = if self.linear_fallback_to_lm != Some(false) {
                AUTOBKSolver::LegacyLm
            } else {
                default_autobk_fallback_solver()
            };
            self.linear_fallback_solver = Some(fallback_solver);
        }

        if self.linear_fallback_to_lm.is_none() {
            self.linear_fallback_to_lm = Some(true);
        }

        if self.linear_workspace_cache.is_none() {
            self.linear_workspace_cache = Some(true);
        }

        Ok(())
    }

    fn validate_input_vectors(
        energy: &DVector<f64>,
        mu: &DVector<f64>,
    ) -> Result<(), BackgroundError> {
        if energy.len() != mu.len() {
            return Err(DataError::LengthMismatch {
                energy_len: energy.len(),
                mu_len: mu.len(),
            }
            .into());
        }
        if energy.len() < 2 {
            return Err(DataError::InsufficientData {
                min: 2,
                actual: energy.len(),
            }
            .into());
        }

        let non_finite = energy
            .iter()
            .zip(mu.iter())
            .enumerate()
            .filter_map(|(index, (e, m))| (!e.is_finite() || !m.is_finite()).then_some(index))
            .collect::<Vec<_>>();
        if !non_finite.is_empty() {
            return Err(DataError::NonFiniteValues {
                indices: non_finite,
            }
            .into());
        }

        for index in 1..energy.len() {
            let prev = energy[index - 1];
            let curr = energy[index];
            if curr < prev {
                return Err(DataError::NonMonotonicEnergy { index, prev, curr }.into());
            }
        }

        Ok(())
    }

    fn build_knot_domain(spl_k: &DVector<f64>, order: usize) -> Result<Vec<f64>, BackgroundError> {
        if spl_k.is_empty() {
            return Err(BackgroundError::SplineKnotsFailed {
                kmin: 0.0,
                kmax: 0.0,
            });
        }

        let nspl = spl_k.len();
        if nspl <= order + 1 {
            return Err(BackgroundError::Other {
                message: format!(
                    "insufficient spline control points for order {}: {}",
                    order, nspl
                ),
            });
        }

        let qmin = spl_k[0];
        let qmax = spl_k[nspl - 1];
        if !qmin.is_finite() || !qmax.is_finite() || qmax <= qmin {
            return Err(BackgroundError::SplineKnotsFailed {
                kmin: qmin,
                kmax: qmax,
            });
        }

        let mut knots = Vec::with_capacity(nspl + order);
        for i in 0..order {
            knots.push(qmin - 1e-4 * (order - i) as f64);
        }
        let interior_count = nspl - order;
        if interior_count < 2 {
            return Err(BackgroundError::Other {
                message: format!(
                    "insufficient interior knots for order {} and control points {}",
                    order, nspl
                ),
            });
        }
        for j in 0..interior_count {
            knots.push(qmin + j as f64 * (qmax - qmin) / (interior_count - 1) as f64);
        }
        let qlast = *knots.last().ok_or_else(|| BackgroundError::Other {
            message: "failed to build interior knot domain".to_string(),
        })?;
        for i in 0..order {
            knots.push(qlast + 1e-4 * (i + 1) as f64);
        }

        Ok(knots)
    }

    fn solve_lm_problem(problem: AUTOBKSpline) -> Result<AUTOBKSpline, BackgroundError> {
        let (fit_result, report) = LevenbergMarquardt::new()
            .with_gtol(1.0e-5)
            .with_ftol(1.0e-5)
            .with_xtol(1.0e-5)
            .with_stepbound(1.0e-5)
            .minimize(problem);
        if !report.termination.was_successful() || !report.objective_function.is_finite() {
            return Err(BackgroundError::OptimizationFailed {
                reason: format!(
                    "legacy LM terminated with {:?} after {} evaluations",
                    report.termination, report.number_of_evaluations
                ),
            });
        }
        Ok(fit_result)
    }

    #[cfg(feature = "trust-region")]
    fn solve_trust_region_problem(problem: AUTOBKSpline) -> Result<AUTOBKSpline, BackgroundError> {
        let residual_len = problem.residual_vec(&problem.coefs).len();
        let factor = AUTOBKSplineFactor {
            spline: problem.clone(),
            residual_len,
        };

        let mut apex_problem = ApexProblem::new(ApexJacobianMode::Dense);
        let coefs_key = apex_problem.add_variable(
            ManifoldType::RN,
            ApexDVector::from_column_slice(problem.coefs.as_slice()),
        );
        apex_problem.add_residual_block(&[coefs_key], Box::new(factor), None);

        let config = DogLegConfig::new()
            .with_linear_solver_type(LinearSolverType::DenseQR)
            .with_max_iterations(100)
            .with_cost_tolerance(1.0e-10)
            .with_parameter_tolerance(1.0e-10)
            .with_gradient_tolerance(1.0e-10);
        let mut optimizer = DogLeg::with_config(config);
        let result = optimizer.optimize(&mut apex_problem).map_err(|err| {
            BackgroundError::OptimizationFailed {
                reason: format!("trust-region DogLeg solver failed: {err}"),
            }
        })?;

        let coefs = result
            .parameters
            .get(coefs_key)
            .map(|variable| variable.to_dvector())
            .ok_or_else(|| BackgroundError::OptimizationFailed {
                reason: "trust-region result missing spline coefficients".to_string(),
            })?;

        if coefs.len() != problem.coefs.len() {
            return Err(BackgroundError::OptimizationFailed {
                reason: format!(
                    "trust-region returned {} coefficients, expected {}",
                    coefs.len(),
                    problem.coefs.len()
                ),
            });
        }

        let mut fit_result = problem;
        fit_result.coefs = DVector::from_vec(coefs.as_slice().to_vec());
        Ok(fit_result)
    }

    #[cfg(not(feature = "trust-region"))]
    fn solve_trust_region_problem(_problem: AUTOBKSpline) -> Result<AUTOBKSpline, BackgroundError> {
        Err(BackgroundError::OptimizationFailed {
            reason: "trust-region DogLeg solver requires the trust-region feature".to_string(),
        })
    }

    fn solve_nonlinear_problem(
        problem: AUTOBKSpline,
        solver: AUTOBKSolver,
    ) -> Result<AUTOBKSpline, BackgroundError> {
        match solver {
            AUTOBKSolver::TrustRegionDogLeg => Self::solve_trust_region_problem(problem),
            AUTOBKSolver::LegacyLm => Self::solve_lm_problem(problem),
            AUTOBKSolver::LinearDirect => Err(BackgroundError::DirectSolverFailed {
                reason: "LinearDirect cannot be used as the fallback solver for LinearDirect"
                    .to_string(),
            }),
        }
    }

    fn solve_spline_problem(
        &self,
        spline_opt: AUTOBKSpline,
    ) -> Result<AUTOBKSpline, BackgroundError> {
        match self.solver.unwrap_or(AUTOBKSolver::LinearDirect) {
            AUTOBKSolver::TrustRegionDogLeg => Self::solve_trust_region_problem(spline_opt),
            AUTOBKSolver::LegacyLm => Self::solve_lm_problem(spline_opt),
            AUTOBKSolver::LinearDirect => {
                let clamp_policy = self
                    .clamp_scale_policy
                    .unwrap_or(AUTOBKClampScalePolicy::Fixed);
                let regularization = self
                    .linear_regularization
                    .unwrap_or(DEFAULT_LINEAR_REGULARIZATION)
                    .max(0.0);
                let condition_limit = self
                    .linear_condition_limit
                    .unwrap_or(DEFAULT_LINEAR_CONDITION_LIMIT)
                    .max(1.0);
                let residual_ratio_limit = self
                    .linear_residual_ratio_limit
                    .unwrap_or(DEFAULT_LINEAR_RESIDUAL_RATIO_LIMIT)
                    .max(1.0);
                let use_workspace_cache = self.linear_workspace_cache.unwrap_or(true);
                let fallback_enabled = self.linear_fallback_to_lm.unwrap_or(true);
                let fallback_solver = self
                    .linear_fallback_solver
                    .unwrap_or_else(default_autobk_fallback_solver);

                match spline_opt.solve_linear_direct(
                    clamp_policy,
                    regularization,
                    condition_limit,
                    residual_ratio_limit,
                    use_workspace_cache,
                ) {
                    Ok(coefs) => {
                        let mut fit_result = spline_opt;
                        fit_result.coefs = coefs;
                        Ok(fit_result)
                    }
                    Err(linear_err) => {
                        if fallback_enabled {
                            // LinearDirect skips building the basis, so the
                            // fallback would otherwise recompute it on every
                            // iteration. Build it once here instead.
                            let mut spline_opt = spline_opt;
                            spline_opt.ensure_precomputed_basis();
                            Self::solve_nonlinear_problem(spline_opt, fallback_solver)
                        } else {
                            Err(linear_err)
                        }
                    }
                }
            }
        }
    }

    /// Fit a smooth background from borrowed energy (eV) and absorption arrays.
    ///
    /// Calculates missing normalization, mutates the supplied normalization settings
    /// and replaces this object's cached outputs: bkg in absorption units and
    /// dimensionless chie on the energy grid, plus dimensionless unweighted chi on
    /// the zero-origin k grid in Å⁻¹. Returns this object for chaining. Scalar
    /// defaults and the resolved ek0 are retained; automatic kmax/nknots remain
    /// local. Earlier parameter filling and normalization can remain after failure.
    ///
    /// Inputs must be finite, matched and nondecreasing with enough post-edge
    /// points; duplicate energies may be nudged in an internal copy. Invalid
    /// cutoffs, spline geometry, unavailable solver features or failed solves
    /// return `BackgroundError`. FixedPenalty requires LinearDirect and never
    /// switches to a legacy fallback objective. Input arrays are unchanged.
    pub fn calc_background(
        &mut self,
        energy: &DVector<f64>,
        mu: &DVector<f64>,
        normalization_param: &mut Option<normalization::NormalizationMethod>,
    ) -> Result<&mut Self, BackgroundError> {
        // Fill in default values for parameters that are not set
        self.fill_parameter()?;
        Self::validate_input_vectors(energy, mu)?;
        let fixed_penalty = self.clamp_scale_policy == Some(AUTOBKClampScalePolicy::FixedPenalty);
        if fixed_penalty {
            if self.solver != Some(AUTOBKSolver::LinearDirect) {
                return Err(BackgroundError::DirectSolverFailed {
                    reason: "FixedPenalty requires LinearDirect; select Fixed or TwoPass for legacy solvers".into(),
                });
            }
            fixed::validate(self)?;
        }

        let rbkg = self.rbkg.ok_or_else(|| DataError::MissingData {
            field: "rbkg".to_string(),
        })?;
        if !rbkg.is_finite() || rbkg <= 0.0 {
            return Err(BackgroundError::InvalidRbkg { rbkg });
        }

        let kmin = self.kmin.ok_or_else(|| DataError::MissingData {
            field: "kmin".to_string(),
        })?;
        let kstep = self.kstep.ok_or_else(|| DataError::MissingData {
            field: "kstep".to_string(),
        })?;
        if !kstep.is_finite() || kstep <= 0.0 {
            return Err(BackgroundError::Other {
                message: format!("invalid kstep: {}", kstep),
            });
        }
        let nfft = self.nfft.ok_or_else(|| DataError::MissingData {
            field: "nfft".to_string(),
        })?;
        if nfft <= 0 {
            return Err(BackgroundError::Other {
                message: format!("invalid nfft: {}", nfft),
            });
        }
        let kweight = self.kweight.ok_or_else(|| DataError::MissingData {
            field: "kweight".to_string(),
        })?;
        let nclamp = self.nclamp.ok_or_else(|| DataError::MissingData {
            field: "nclamp".to_string(),
        })?;
        let clamp_lo = self.clamp_lo.ok_or_else(|| DataError::MissingData {
            field: "clamp_lo".to_string(),
        })?;
        let clamp_hi = self.clamp_hi.ok_or_else(|| DataError::MissingData {
            field: "clamp_hi".to_string(),
        })?;
        let energy_dv = energy;
        let mu_dv = mu;

        let energy = energy_dv;
        let mu = mu_dv;

        // Perform normalization if necessary

        if normalization_param.is_none() {
            let mut normalization_method = normalization::PrePostEdge::new();
            normalization_method.set_e0(self.ek0);
            *normalization_param = Some(normalization::NormalizationMethod::PrePostEdge(
                normalization_method,
            ));
        }
        let normalization_method =
            normalization_param
                .as_mut()
                .ok_or_else(|| DataError::MissingData {
                    field: "normalization method".to_string(),
                })?;

        if let Some(ek0) = self.ek0 {
            if ek0 < energy_dv.min() || ek0 > energy_dv.max() {
                self.ek0 = None;
                normalization_method.set_e0(None);
            }
        }

        let e0 = normalization_method.get_e0();
        let mut edge_step = normalization_method.get_edge_step();
        let ek0 = self.ek0;

        if (ek0.is_none() && e0.is_none()) || edge_step.is_none() {
            normalization_method.normalize(energy_dv, mu_dv)?;
            edge_step = normalization_method.get_edge_step();
        }

        let ek0 = if self.ek0.is_none() {
            normalization_method.get_e0()
        } else {
            self.ek0
        }
        .ok_or_else(|| DataError::MissingData {
            field: "ek0/e0".to_string(),
        })?;
        self.ek0 = Some(ek0);

        let edge_step = edge_step
            .or_else(|| normalization_method.get_edge_step())
            .ok_or_else(|| DataError::MissingData {
                field: "edge_step".to_string(),
            })?;
        if !edge_step.is_finite() {
            return Err(BackgroundError::Other {
                message: format!("edge_step is not finite: {}", edge_step),
            });
        }
        if fixed_penalty && edge_step <= 0.0 {
            return Err(BackgroundError::Other {
                message: "FixedPenalty requires a positive edge step".into(),
            });
        }
        let edge_step = if fixed_penalty {
            edge_step
        } else {
            edge_step.max(1.0e-12)
        };
        // Rbkg Algorithm
        let energy_slice = energy.as_slice();
        let iek0 = mathutils::index_of_sorted(energy_slice, &ek0)?;
        let mut rgrid = std::f64::consts::PI / (kstep * nfft as f64);

        if rbkg < (2.0 * rgrid) {
            rgrid *= 2.0;
        }

        let enpe = DVector::from_iterator(
            energy.len() - iek0,
            energy.iter().skip(iek0).map(|x| x - ek0),
        );
        let kraw = enpe.map(|x| x.signum() * (xafsutils::constants::ETOK * x.abs()).sqrt());

        let kmax = self
            .kmax
            .map(|max_value| max_value.min(kraw.max()).max(0.0))
            .unwrap_or_else(|| kraw.max());
        if kmax <= kmin {
            return Err(BackgroundError::SplineKnotsFailed { kmin, kmax });
        }

        let kout = dvector_arange(0.0, (1.01 + kmax / kstep).floor(), 1.0) * kstep;
        if kout.len() < 2 {
            return Err(DataError::InsufficientData {
                min: 2,
                actual: kout.len(),
            }
            .into());
        }

        let iemax = energy.len().min(
            2 + mathutils::index_of_sorted(
                energy_slice,
                &(ek0 + kmax.powi(2) / xafsutils::constants::ETOK),
            )?,
        ) - 1;
        if iemax <= iek0 {
            return Err(DataError::InsufficientData {
                min: 3,
                actual: iemax.saturating_sub(iek0) + 1,
            }
            .into());
        }

        let chi_std = match (&self.k_std, &self.chi_std) {
            (Some(k_std), Some(chi_std)) => {
                Some(kout.interpolate(k_std.as_slice(), chi_std.as_slice())?)
            }
            (None, None) => None,
            _ => {
                return Err(DataError::MissingData {
                    field: "k_std and chi_std must both be set".to_string(),
                }
                .into())
            }
        };

        let ftwin = kout
            .map(|x| x.powi(kweight))
            .component_mul(&xafsutils::ftwindow(
                &kout,
                Some(kmin),
                Some(kmax),
                self.dk,
                self.dk,
                Some(self.window),
            )?);

        let mut nspl = 1 + (2.0 * rbkg * (kmax - kmin) / std::f64::consts::PI).floor() as i32;
        let cutoff = 1.0 + (nspl - 1) as f64 * std::f64::consts::PI / (2.0 * rgrid * (kmax - kmin));
        let irbkg = if fixed_penalty {
            cutoff.floor()
        } else {
            cutoff.round()
        } as i32;

        if let Some(nknots) = self.nknots {
            nspl = nknots;
        }

        nspl = nspl.clamp(5, 128);

        // !todo!("Finish implementing this part of the code");
        let mut spl_y = DVector::from_element(nspl as usize, 1.0);
        let mut spl_k = DVector::zeros(nspl as usize);
        let kraw_slice = kraw.as_slice();
        for i in 0..nspl as usize {
            let q = kmin + i as f64 * (kmax - kmin) / (nspl - 1) as f64;
            let ik = mathutils::index_nearest_sorted(kraw_slice, &q)?;
            let i1 = (ik + 5).min(kraw.len() - 1);
            let i2 = (ik as i32 - 5).max(0) as usize;
            spl_k[i] = kraw[ik];
            spl_y[i] = (2.0 * mu[ik + iek0] + mu[i1 + iek0] + mu[i2 + iek0]) / 4.0;
        }

        let order = 3;
        // Validate knot-domain construction used for AUTOBK bounds reasoning.
        let _knot_domain = Self::build_knot_domain(&spl_k, order)?;
        let (knots, coefs) = if fixed_penalty {
            let knots = spline::cubic_knots(spl_k.as_slice())
                .map_err(|message| BackgroundError::Other { message })?;
            (knots, vec![0.0; nspl as usize])
        } else {
            spline::interpolate(spl_k.as_slice(), spl_y.as_slice(), order)
                .map_err(|message| BackgroundError::Other { message })?
        };

        // Calculate the mu interpolated to the k grid
        let kraw_fit = kraw.rows(0, iemax - iek0 + 1).into_owned();
        let mu_fit = mu.rows(iek0, iemax - iek0 + 1).into_owned();
        let mu_out = if fixed_penalty {
            DVector::from_vec(
                spline::cubic_resample(kraw_fit.as_slice(), mu_fit.as_slice(), kout.as_slice())
                    .map_err(|message| BackgroundError::Other { message })?,
            ) / edge_step
        } else {
            kout.interpolate(kraw_fit.as_slice(), mu_fit.as_slice())?
        };

        // Built once here rather than on every LM iteration: the basis depends
        // only on knots, order and evaluation points, all fixed for the
        // lifetime of this spline. See `residual_jacobian_with_scale`.
        //
        // Only the iterative solvers consume it. LinearDirect solves through
        // `linear_design_matrix` and never asks for this Jacobian, so building
        // it there is pure overhead — measured at +20% on that path.
        let precomputed_basis = match self.solver.unwrap_or(AUTOBKSolver::LinearDirect) {
            AUTOBKSolver::LegacyLm | AUTOBKSolver::TrustRegionDogLeg => -splev_jacobian(
                knots.clone(),
                coefs.clone(),
                order,
                kout.data.as_vec().clone(),
                3,
            ),
            AUTOBKSolver::LinearDirect => DMatrix::zeros(0, 0),
        };

        let spline_opt = AUTOBKSpline {
            coefs: DVector::from_vec(coefs),
            knots: DVector::from_vec(knots),
            order,
            precomputed_basis,
            irbkg: irbkg.max(1) as usize,
            nfft: nfft as usize,
            kraw: kraw_fit,
            mu: mu_out,
            kout: kout.clone(),
            ftwin,
            kweight,
            chi_std,
            nclamp,
            clamp_lo,
            clamp_hi,
            kstep,
            ..Default::default()
        };

        if fixed_penalty {
            let (bkg, chi) = fixed::solve(&spline_opt, self)?;
            let mut obkg = mu.clone();
            for (i, value) in bkg.iter().enumerate() {
                obkg[iek0 + i] = value * edge_step;
            }
            self.chie = Some((mu - &obkg) / edge_step);
            self.bkg = Some(obkg);
            self.k = Some(kout);
            self.chi = Some(chi);
            return Ok(self);
        }

        let fit_result = self.solve_spline_problem(spline_opt)?;

        let (bkg, chi) = spline_eval_nalgebra(
            &fit_result.kraw,
            &fit_result.mu,
            &fit_result.knots,
            &fit_result.coefs,
            fit_result.order,
            &fit_result.kout,
        );

        let mut obkg = mu.clone();
        for i in 0..bkg.len() {
            obkg[iek0 + i] = bkg[i];
        }

        self.chie = Some((mu - &obkg) / edge_step);
        self.bkg = Some(obkg);
        self.k = Some(kout);
        self.chi = Some(chi / edge_step);

        Ok(self)
    }

    /// Read the stored value without recalculating.
    ///
    /// Background energy origin in eV. `None` uses the normalization E0. An explicit
    /// value changes conversion to k without replacing a valid normalization E0.
    pub fn get_ek0(&self) -> Option<&f64> {
        self.ek0.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Low-R background cutoff in Å; recommended default 1.0, finite and positive.
    /// Larger values increase spline flexibility and may remove structural signal.
    /// This is an algorithmic cutoff, not an estimate of a bond length.
    pub fn get_rbkg(&self) -> Option<&f64> {
        self.rbkg.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Requested number of fitted cubic-spline coefficients. `None` uses
    /// `1 + floor(2*rbkg*(kmax-kmin)/pi)`; the final count is bounded to 5–128.
    /// Changing it overrides flexibility but does not independently change the low-R cutoff.
    pub fn get_nknots(&self) -> Option<&i32> {
        self.nknots.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Lower spline/window boundary in Å⁻¹; default 0.0. The returned k grid still
    /// starts at zero; endpoint penalties use that output grid, not this boundary.
    pub fn get_kmin(&self) -> Option<&f64> {
        self.kmin.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Upper spline/window boundary in Å⁻¹. `None` uses the available post-edge range;
    /// an explicit value is capped at that range. Must exceed `kmin`.
    pub fn get_kmax(&self) -> Option<&f64> {
        self.kmax.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Uniform output k spacing in Å⁻¹; default 0.05, finite and positive.
    /// This sets the physical R grid. FixedPenalty uses the fixed numerical
    /// amplitude reference `0.05 / sqrt(pi)`; legacy policies instead use the
    /// actual `kstep / sqrt(pi)`.
    pub fn get_kstep(&self) -> Option<&f64> {
        self.kstep.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Number of output k samples penalized at each enabled end; default 3.
    /// FixedPenalty uses up to the available point count and includes the last sample.
    /// Zero disables endpoint rows; FixedPenalty rejects negative counts.
    pub fn get_nclamp(&self) -> Option<&i32> {
        self.nclamp.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Low-k endpoint weight; default 0 (disabled). FixedPenalty uses its absolute
    /// value and squares that weight in the objective, so doubling it quadruples
    /// that end's contribution.
    pub fn get_clamp_lo(&self) -> Option<&i32> {
        self.clamp_lo.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// High-k endpoint weight; default 1. FixedPenalty uses its absolute value;
    /// zero disables that end. The returned chi is not multiplied by this weight.
    pub fn get_clamp_hi(&self) -> Option<&i32> {
        self.clamp_hi.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Internal Fourier length; default 2048. Use enough points for the complete
    /// output grid; the low-level FFT pads short arrays and truncates longer arrays.
    pub fn get_nfft(&self) -> Option<&i32> {
        self.nfft.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Objective-only standard, interpolated from `k_std` onto the output k grid.
    /// FixedPenalty expects dimensionless chi. Legacy policies subtract the
    /// standard before edge-step division, so they require values in original
    /// mu units. It is never subtracted from the returned chi.
    pub fn get_chi_std(&self) -> Option<&DVector<f64>> {
        self.chi_std.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Wave-number grid of `chi_std`, in Å⁻¹; both arrays must be supplied together.
    /// Use finite matching arrays with strictly increasing k for interpolation.
    pub fn get_k_std(&self) -> Option<&DVector<f64>> {
        self.k_std.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Integer k-weight used in the background objective; default 1.
    /// It is separate from public forward-transform weighting and does not multiply
    /// the returned unweighted chi. Larger weights emphasize high-k data and noise.
    pub fn get_kweight(&self) -> Option<&i32> {
        self.kweight.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Window family for the background objective; default Hanning.
    /// Changing the family alters the objective and can affect the fitted background.
    pub fn get_window(&self) -> FTWindow {
        self.window
    }

    /// Read the stored value without recalculating.
    ///
    /// Internal window parameter; default 0.1, applied at both ends.
    /// Usually a taper width in Å⁻¹; Kaiser–Bessel also uses its value for shape,
    /// Gaussian uses it for width, and fractional Hanning uses fractional geometry.
    pub fn get_dk(&self) -> Option<&f64> {
        self.dk.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Smooth background on the original energy grid, in input mu units.
    /// Outside the fitted interval this equals the original mu; `None` before calculation.
    pub fn get_bkg(&self) -> Option<&DVector<f64>> {
        self.bkg.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Dimensionless `(mu - bkg) / edge_step` on the original energy grid.
    /// It is zero outside the fitted interval; `None` before calculation.
    pub fn get_chie(&self) -> Option<&DVector<f64>> {
        self.chie.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Uniform output wave-number grid in Å⁻¹ starting at zero.
    /// `None` before calculation.
    pub fn get_k(&self) -> Option<&DVector<f64>> {
        self.k.as_ref()
    }

    /// Read the stored value without recalculating.
    ///
    /// Dimensionless unweighted EXAFS on `k`, without objective-only standard
    /// subtraction. `None` before calculation.
    pub fn get_chi(&self) -> Option<&DVector<f64>> {
        self.chi.as_ref()
    }

    /// Allocate chi multiplied by the background k-weight; no window is applied.
    /// Units are Å⁻ʷ for dimensionless chi and k-weight w; `None` before calculation.
    pub fn get_chi_kweighted(&self) -> Option<DVector<f64>> {
        let kweight = self.kweight?;
        let k = self.k.clone()?;
        let chi = self.chi.clone()?;

        if kweight == 0 {
            Some(chi)
        } else {
            Some(chi.component_mul(&k.map(|x| x.powi(kweight))))
        }
    }

    /// Calculate a dimensionless background window on the stored k grid.
    /// No k-weight is included. Returns `None` without a grid or if window creation fails.
    pub fn get_ftwin(&self) -> Option<DVector<f64>> {
        let k = self.k.as_ref()?;

        let ftwin =
            xafsutils::ftwindow(k, self.kmin, self.kmax, self.dk, self.dk, Some(self.window))
                .ok()?;

        Some(ftwin)
    }
}

/// Evaluation of the spline used in AUTOBK
///
/// Inputs and outputs use nalgebra vectors.
///
/// # Arguments
///
/// * `kraw` - kraw, the k grid converted from energy
/// * `mu` - mu(E)
/// * `knots` - knots of the spline
/// * `coefs` - coefficients of the spline
/// * `order` - order of the spline
/// * `kout` - k grid ready for FFT
fn spline_eval_nalgebra(
    kraw: &DVector<f64>,
    mu: &DVector<f64>,
    knots: &DVector<f64>,
    coefs: &DVector<f64>,
    order: usize,
    kout: &DVector<f64>,
) -> (DVector<f64>, DVector<f64>) {
    let knots_vec = knots.as_slice().to_vec();
    let coefs_vec = coefs.as_slice().to_vec();

    let bkg = DVector::from_vec(spline::evaluate(
        &knots_vec,
        &coefs_vec,
        order,
        kraw.as_slice(),
        true,
    ));

    // experimental
    let bkg_out = DVector::from_vec(spline::evaluate(
        &knots_vec,
        &coefs_vec,
        order,
        kout.as_slice(),
        true,
    ));

    let chi = mu - &bkg_out;

    (bkg, chi.clone())
}

/// Struct for solving Levenberg-Marquardt optimization for AUTOBK
#[derive(Debug, Clone, PartialEq)]
struct AUTOBKSpline {
    pub coefs: DVector<f64>,
    pub knots: DVector<f64>,
    pub order: usize,
    pub irbkg: usize,
    pub nfft: usize,
    pub kraw: DVector<f64>,
    pub mu: DVector<f64>,
    pub kout: DVector<f64>,
    pub ftwin: DVector<f64>,
    pub kweight: i32,
    pub chi_std: Option<DVector<f64>>,
    pub nclamp: i32,
    pub clamp_lo: i32,
    pub clamp_hi: i32,
    pub kstep: f64,
    pub scale: f64,
    /// B-spline basis Jacobian (`kout.len()` x `coefs.len()`), negated.
    /// Depends only on knots, order and evaluation points, so it is computed
    /// once instead of on every Levenberg-Marquardt iteration. Empty when the
    /// spline was default-constructed; see `residual_jacobian_with_scale`.
    pub precomputed_basis: DMatrix<f64>,
}

impl Default for AUTOBKSpline {
    fn default() -> Self {
        AUTOBKSpline {
            coefs: DVector::zeros(0),
            knots: DVector::zeros(0),
            order: 3,
            irbkg: 1,
            nfft: 2048,
            kraw: DVector::zeros(0),
            mu: DVector::zeros(0),
            kout: DVector::zeros(0),
            ftwin: DVector::zeros(0),
            kweight: 1,
            chi_std: None,
            nclamp: 0,
            clamp_lo: 1,
            clamp_hi: 1,
            kstep: 0.05,
            scale: 1.0,
            precomputed_basis: DMatrix::zeros(0, 0),
        }
    }
}

impl AUTOBKSpline {
    /// Whether `precomputed_basis` matches this spline's current shape.
    /// A `Default`-constructed spline has an empty basis and fails this.
    fn basis_cache_is_valid(&self) -> bool {
        self.precomputed_basis.nrows() == self.kout.len()
            && self.precomputed_basis.ncols() == self.coefs.len()
    }

    /// Build `precomputed_basis` unless it is already valid. Used when a
    /// solver that skipped the precomputation hands the spline to one that
    /// needs it.
    fn ensure_precomputed_basis(&mut self) {
        if self.basis_cache_is_valid() {
            return;
        }
        self.precomputed_basis = -splev_jacobian(
            self.knots.data.as_vec().clone(),
            self.coefs.data.as_vec().clone(),
            self.order,
            self.kout.data.as_vec().clone(),
            3,
        );
    }

    fn chi_for_coefs(&self, coefs: &DVector<f64>) -> DVector<f64> {
        let (_, chi) = spline_eval_nalgebra(
            &self.kraw,
            &self.mu,
            &self.knots,
            coefs,
            self.order,
            &self.kout,
        );

        if let Some(chi_std) = self.chi_std.as_ref() {
            chi - chi_std
        } else {
            chi
        }
    }

    fn fft_residual_head(&self, chi: &DVector<f64>) -> DVector<f64> {
        let fft = chi
            .component_mul(&self.ftwin)
            .xftf_fast(self.nfft, self.kstep);
        let upper = self.irbkg.min(fft.len());
        fft[..upper].realimg()
    }

    fn clamp_len(&self, vec_len: usize) -> usize {
        if self.nclamp <= 0 {
            return 0;
        }
        (self.nclamp as usize).min(vec_len.saturating_sub(1))
    }

    fn clamp_scale(&self, coefs: &DVector<f64>) -> f64 {
        if self.nclamp == 0 {
            return 1.0;
        }
        let chi = self.chi_for_coefs(coefs);
        let out = self.fft_residual_head(&chi);
        if out.is_empty() {
            return 1.0;
        }

        1.0 + 100.0 * out.dot(&out) / out.len() as f64
    }

    /// Build the historical residual with dynamic or explicitly frozen clamp scale.
    /// The head contains real/imaginary FFT entries scaled by actual kstep/sqrt(pi).
    /// Dynamic s(c)=1+100*mean(head(c)^2); endpoint rows are weight*s(c)*chi(c).
    /// Both ends contribute rows, including zero-weight rows, and the high-end
    /// slice excludes the final chi sample. This is not FixedPenalty.
    pub fn residual_vec_with_scale(
        &self,
        coefs: &DVector<f64>,
        clamp_scale_override: Option<f64>,
    ) -> DVector<f64> {
        let chi = self.chi_for_coefs(coefs);
        let mut out = self.fft_residual_head(&chi);

        if self.nclamp == 0 {
            return out;
        }

        let nclamp = self.clamp_len(chi.len());
        if nclamp == 0 {
            return out;
        }

        let scale = clamp_scale_override.unwrap_or_else(|| self.clamp_scale(coefs));
        let low_clamp = self.clamp_lo as f64 * scale * chi.view((0, 0), (nclamp, 1));
        let high_start = chi.len() - nclamp - 1;
        let high_clamp = self.clamp_hi as f64 * scale * chi.view((high_start, 0), (nclamp, 1));

        out.extend(low_clamp.data.as_vec().to_owned());
        out.extend(high_clamp.data.as_vec().to_owned());
        out
    }

    pub fn residual_vec(&self, coefs: &DVector<f64>) -> DVector<f64> {
        self.residual_vec_with_scale(coefs, None)
    }

    /// Differentiate the historical residual with respect to spline coefficients.
    /// For dynamic clamps the product rule gives
    /// d(weight*s*chi)/dc = weight*(s*dchi/dc + chi*ds/dc), with
    /// ds/dc = 200*dot(head,dhead/dc)/len(head). A frozen scale has ds/dc=0.
    /// The complete term matches the actual residual since 0.2.4.
    pub fn residual_jacobian_with_scale(
        &self,
        coefs: &DVector<f64>,
        clamp_scale_override: Option<f64>,
    ) -> DMatrix<f64> {
        // For a dynamic clamp s(c) * chi(c), include both terms of the
        // product rule. A frozen scale is used by the linear solver and has
        // zero derivative. Compute the shared residual only once per Jacobian.
        let dynamic =
            (self.clamp_len(self.kout.len()) > 0 && clamp_scale_override.is_none()).then(|| {
                let chi = self.chi_for_coefs(coefs);
                let head = self.fft_residual_head(&chi);
                (chi, head)
            });
        let scale = clamp_scale_override.unwrap_or_else(|| {
            dynamic.as_ref().map_or(1.0, |(_, head)| {
                if head.is_empty() {
                    1.0
                } else {
                    1.0 + 100.0 * head.norm_squared() / head.len() as f64
                }
            })
        });

        // `splev_jacobian` reads `coefs` only for its length (it sizes the
        // output and bounds-checks the write); every entry comes from
        // the B-spline basis over the knots, order and evaluation points. Those three
        // are fixed for the lifetime of the spline — `set_params` writes only
        // `coefs` — so this matrix is identical on every LM iteration and is
        // built once in `calc_background`.
        //
        // The guard keeps a `Default`-constructed spline (which has an empty
        // basis) correct rather than silently wrong.
        let fallback;
        let spline_jacobian: &DMatrix<f64> = if self.basis_cache_is_valid() {
            &self.precomputed_basis
        } else {
            fallback = -splev_jacobian(
                self.knots.data.as_vec().clone(),
                self.coefs.data.as_vec().clone(),
                self.order,
                self.kout.data.as_vec().clone(),
                3,
            );
            &fallback
        };

        let jacobian_columns = spline_jacobian
            .column_iter()
            .map(|chi_der| {
                let fft = chi_der
                    .component_mul(&self.ftwin)
                    .xftf_fast(self.nfft, self.kstep);
                let upper = self.irbkg.min(fft.len());
                let mut out: DVector<f64> = fft[..upper].realimg();

                if self.nclamp == 0 {
                    return out;
                }

                let nclamp = self.clamp_len(chi_der.len());
                if nclamp == 0 {
                    return out;
                }

                let scale_der = dynamic.as_ref().map_or(0.0, |(_, head)| {
                    if head.is_empty() {
                        0.0
                    } else {
                        200.0 * head.dot(&out) / head.len() as f64
                    }
                });
                let high_start = chi_der.len() - nclamp - 1;
                for (start, weight) in [(0, self.clamp_lo), (high_start, self.clamp_hi)] {
                    out.extend((start..start + nclamp).map(|i| {
                        let chi = dynamic.as_ref().map_or(0.0, |(chi, _)| chi[i]);
                        weight as f64 * (scale * chi_der[i] + scale_der * chi)
                    }));
                }
                out
            })
            .collect::<Vec<DVector<f64>>>();

        DMatrix::from_columns(&jacobian_columns)
    }

    pub fn residual_jacobian(&self, coefs: &DVector<f64>) -> DMatrix<f64> {
        self.residual_jacobian_with_scale(coefs, None)
    }

    fn workspace_signature(&self, clamp_scale: f64) -> u64 {
        let mut hasher = DefaultHasher::new();
        hasher.write_u64(self.order as u64);
        hasher.write_u64(self.irbkg as u64);
        hasher.write_u64(self.nfft as u64);
        hasher.write_u64(self.nclamp.max(0) as u64);
        hasher.write_u64(self.clamp_lo.max(0) as u64);
        hasher.write_u64(self.clamp_hi.max(0) as u64);
        hasher.write_u64(self.kweight.max(0) as u64);
        hasher.write_u64(clamp_scale.to_bits());
        hasher.write_u64(self.kstep.to_bits());
        hasher.write_u64(self.knots.len() as u64);
        hasher.write_u64(self.coefs.len() as u64);
        hasher.write_u64(self.kout.len() as u64);
        hasher.write_u64(self.ftwin.len() as u64);

        for value in self.knots.iter() {
            hasher.write_u64(value.to_bits());
        }
        for value in self.kout.iter() {
            hasher.write_u64(value.to_bits());
        }
        for value in self.ftwin.iter() {
            hasher.write_u64(value.to_bits());
        }

        hasher.finish()
    }

    fn linear_design_matrix(
        &self,
        clamp_scale: f64,
        use_workspace_cache: bool,
    ) -> Arc<DMatrix<f64>> {
        let signature = self.workspace_signature(clamp_scale);

        if use_workspace_cache {
            if let Ok(mut cache) = autobk_workspace_cache().lock() {
                if let Some(index) = cache.iter().position(|entry| entry.signature == signature) {
                    if let Some(entry) = cache.remove(index) {
                        let matrix = entry.design_matrix.clone();
                        cache.push_front(entry);
                        return matrix;
                    }
                }
            }
        }

        let matrix = Arc::new(self.residual_jacobian_with_scale(&self.coefs, Some(clamp_scale)));
        if use_workspace_cache {
            if let Ok(mut cache) = autobk_workspace_cache().lock() {
                cache.retain(|entry| entry.signature != signature);
                cache.push_front(AutobkLinearWorkspace {
                    signature,
                    design_matrix: matrix.clone(),
                });
                cache.truncate(AUTOBK_WORKSPACE_CACHE_CAPACITY);
            }
        }

        matrix
    }

    fn solve_linear_direct_for_scale(
        &self,
        clamp_scale: f64,
        regularization: f64,
        condition_limit: f64,
        residual_ratio_limit: f64,
        use_workspace_cache: bool,
    ) -> Result<DVector<f64>, BackgroundError> {
        if self.coefs.is_empty() {
            return Err(BackgroundError::DirectSolverFailed {
                reason: "empty coefficient vector".to_string(),
            });
        }

        let base_residual =
            self.residual_vec_with_scale(&DVector::zeros(self.coefs.len()), Some(clamp_scale));
        if base_residual.iter().any(|v| !v.is_finite()) {
            return Err(BackgroundError::DirectSolverFailed {
                reason: "non-finite base residual".to_string(),
            });
        }

        let design = self.linear_design_matrix(clamp_scale, use_workspace_cache);
        if design.iter().any(|v| !v.is_finite()) {
            return Err(BackgroundError::DirectSolverFailed {
                reason: "non-finite design matrix".to_string(),
            });
        }

        let mut regularization_candidates = vec![regularization.max(0.0)];
        for candidate in [1.0e-8, 1.0e-6, 1.0e-4, 1.0e-2, 1.0] {
            if candidate > regularization.max(0.0) {
                regularization_candidates.push(candidate);
            }
        }

        let mut last_error = None;

        for reg in regularization_candidates {
            let n_rows = design.nrows();
            let n_cols = design.ncols();
            let reg_sqrt = reg.sqrt();

            let mut augmented_design = DMatrix::zeros(n_rows + n_cols, n_cols);
            for i in 0..n_rows {
                for j in 0..n_cols {
                    augmented_design[(i, j)] = design[(i, j)];
                }
            }
            for i in 0..n_cols {
                augmented_design[(n_rows + i, i)] = reg_sqrt;
            }

            let mut augmented_rhs = DVector::zeros(n_rows + n_cols);
            for i in 0..n_rows {
                augmented_rhs[i] = -base_residual[i];
            }

            let mut column_scales = DVector::from_element(n_cols, 1.0);
            let mut scaled_design = augmented_design;
            for j in 0..n_cols {
                let norm = scaled_design.column(j).norm();
                let scale = if norm.is_finite() && norm > f64::EPSILON {
                    norm
                } else {
                    1.0
                };
                column_scales[j] = scale;
                for i in 0..scaled_design.nrows() {
                    scaled_design[(i, j)] /= scale;
                }
            }

            let svd = scaled_design.svd(true, true);
            let mut sigma_max = 0.0_f64;
            let mut sigma_min = f64::INFINITY;
            for sigma in svd.singular_values.iter() {
                let value = sigma.abs();
                sigma_max = sigma_max.max(value);
                if value > f64::EPSILON {
                    sigma_min = sigma_min.min(value);
                }
            }

            let condition_proxy = if sigma_min.is_finite() && sigma_min > f64::EPSILON {
                sigma_max / sigma_min
            } else {
                f64::INFINITY
            };

            if !condition_proxy.is_finite() || condition_proxy > condition_limit {
                last_error = Some(BackgroundError::DirectSolverIllConditioned {
                    condition_proxy,
                    limit: condition_limit,
                });
                continue;
            }

            let solved_scaled = match svd.solve(&augmented_rhs, f64::EPSILON.sqrt()) {
                Ok(solution) => solution,
                Err(reason) => {
                    last_error = Some(BackgroundError::DirectSolverFailed {
                        reason: format!("svd solve failed at regularization {}: {}", reg, reason),
                    });
                    continue;
                }
            };

            let mut coefs = solved_scaled;
            for j in 0..coefs.len() {
                coefs[j] /= column_scales[j];
            }

            if coefs.iter().any(|v| !v.is_finite()) {
                last_error = Some(BackgroundError::DirectSolverFailed {
                    reason: "non-finite coefficient solution".to_string(),
                });
                continue;
            }

            let solved_residual = &base_residual + design.as_ref() * &coefs;
            let base_norm = base_residual.norm();
            let solved_norm = solved_residual.norm();
            if !base_norm.is_finite() || !solved_norm.is_finite() {
                last_error = Some(BackgroundError::DirectSolverFailed {
                    reason: "non-finite residual norms".to_string(),
                });
                continue;
            }

            if base_norm > f64::EPSILON {
                let ratio = solved_norm / base_norm;
                if !ratio.is_finite() || ratio > residual_ratio_limit {
                    last_error = Some(BackgroundError::DirectSolverFailed {
                        reason: format!(
                            "residual quality ratio {} exceeded limit {} at regularization {}",
                            ratio, residual_ratio_limit, reg
                        ),
                    });
                    continue;
                }
            }

            return Ok(coefs);
        }

        Err(
            last_error.unwrap_or_else(|| BackgroundError::DirectSolverFailed {
                reason: "direct solve failed for all regularization attempts".to_string(),
            }),
        )
    }

    pub fn solve_linear_direct(
        &self,
        clamp_policy: AUTOBKClampScalePolicy,
        regularization: f64,
        condition_limit: f64,
        residual_ratio_limit: f64,
        use_workspace_cache: bool,
    ) -> Result<DVector<f64>, BackgroundError> {
        let first_scale = if self.nclamp == 0 {
            1.0
        } else {
            self.clamp_scale(&self.coefs)
        };

        let mut coefs = self.solve_linear_direct_for_scale(
            first_scale,
            regularization,
            condition_limit,
            residual_ratio_limit,
            use_workspace_cache,
        )?;

        if self.nclamp != 0 && clamp_policy == AUTOBKClampScalePolicy::TwoPass {
            let second_scale = self.clamp_scale(&coefs);
            if second_scale.is_finite() && (second_scale - first_scale).abs() > 1.0e-12 {
                coefs = self.solve_linear_direct_for_scale(
                    second_scale,
                    regularization,
                    condition_limit,
                    residual_ratio_limit,
                    use_workspace_cache,
                )?;
            }
        }

        Ok(coefs)
    }
}

#[cfg(feature = "trust-region")]
#[derive(Clone)]
struct AUTOBKSplineFactor {
    spline: AUTOBKSpline,
    residual_len: usize,
}

#[cfg(feature = "trust-region")]
impl ApexFactor for AUTOBKSplineFactor {
    fn linearize(
        &self,
        params: &[&[f64]],
        residual: &mut [f64],
        jacobian: Option<faer::MatMut<'_, f64>>,
    ) {
        let coefs = params
            .first()
            .filter(|param| param.len() == self.spline.coefs.len())
            .map(|param| DVector::from_column_slice(param))
            .unwrap_or_else(|| self.spline.coefs.clone());
        let base = self.spline.residual_vec(&coefs);
        residual.copy_from_slice(base.as_slice());
        if let Some(mut output) = jacobian {
            let values = self.spline.residual_jacobian(&coefs);
            for col in 0..values.ncols() {
                for row in 0..values.nrows() {
                    output[(row, col)] = values[(row, col)];
                }
            }
        }
    }

    fn residual_dim(&self) -> usize {
        self.residual_len
    }
    fn jacobian_shape(&self) -> (usize, usize) {
        (self.residual_len, self.spline.coefs.len())
    }
}

/// Implementation of LeastSquaresProblem trait for AUTOBK algorithm
impl LeastSquaresProblem<f64, Dyn, Dyn> for AUTOBKSpline {
    type ParameterStorage = Owned<f64, Dyn>;
    type ResidualStorage = Owned<f64, Dyn>;
    type JacobianStorage = Owned<f64, Dyn, Dyn>;

    fn set_params(&mut self, coefs: &DVector<f64>) {
        self.coefs.copy_from(coefs);
    }

    fn params(&self) -> DVector<f64> {
        self.coefs.clone()
    }

    fn residuals(&self) -> Option<DVector<f64>> {
        Some(self.residual_vec(&self.coefs))
    }

    /// Jacobian matrix for the Levenberg-Marquardt optimization
    /// Jacobian matrix is calculated by numerical differentiation using foward difference
    fn jacobian(&self) -> Option<DMatrix<f64>> {
        // let residual_vec = |coefs: &DVector<f64>| AUTOBKSpline::residual_vec(&self, &coefs);
        // Some(self.coefs.jacobian(&residual_vec))

        // let start = Instant::now();

        // let jac1 = self.coefs.jacobian(&residual_vec);
        // let duration = start.elapsed();

        // println!("jac1: {}", duration.as_secs_f64());

        // let start = Instant::now();
        // let jac2 = self.residual_jacobian(&self.coefs);

        // let duration = start.elapsed();
        // println!("jac2: {}", duration.as_secs_f64());

        // println!("jac1: {:?}", jac1.shape());
        // println!("jac2: {:?}", jac2.shape());

        // jac1.iter().zip(jac2.iter()).for_each(|(x, y)| {
        //     println!("x: {}, y: {}", x, y);
        //     assert_abs_diff_eq!(x, y, epsilon = 1.0e-1);
        // });

        // Some(jac2)
        Some(self.residual_jacobian(&self.coefs))
    }
}

fn dvector_arange(start: f64, stop: f64, step: f64) -> DVector<f64> {
    if !step.is_finite() || step <= 0.0 {
        return DVector::zeros(0);
    }

    let mut values = Vec::new();
    let mut value = start;
    while value < stop {
        values.push(value);
        value += step;
    }
    DVector::from_vec(values)
}

/// Unimplemented background-removal placeholder.
/// Selecting it returns `BackgroundError::NotImplemented` during calculation.
#[derive(Debug, Clone, PartialEq, Default, Serialize, Deserialize)]
#[serde(default)]
pub struct ILPBkg {}

impl ILPBkg {
    /// Create an empty placeholder, not an implemented background method.
    pub fn new() -> ILPBkg {
        ILPBkg::default()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::xafs::io;
    use crate::xafs::normalization::PrePostEdge;
    use data_reader::reader::{load_txt_f64, Delimiter, ReaderParams};

    use crate::xafs::tests::PARAM_LOADTXT;
    use crate::xafs::tests::TEST_TOL;
    use crate::xafs::tests::TOP_DIR;
    use approx::assert_abs_diff_eq;

    const CHI_MSE_TOL: f64 = 1.0e-4;

    #[test]
    fn test_autobk() -> Result<(), Box<dyn Error>> {
        let acceptable_e0_diff = 1.5;

        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut xafs_test_group = io::load_spectrum_QAS_trans(&path).unwrap();

        // let mut pre_post_edge = PrePostEdge::new();
        // let _ = pre_post_edge.fill_parameter(
        //     &xafs_test_group.energy.clone().unwrap(),
        //     &xafs_test_group.mu.clone().unwrap(),
        // );

        xafs_test_group
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let mut autobk = AUTOBK::new();

        autobk.calc_background(
            &xafs_test_group.energy.clone().unwrap(),
            &xafs_test_group.mu.clone().unwrap(),
            &mut xafs_test_group.normalization,
        )?;

        // Test for chi with larch
        // The chi is not exactly the same as the one calculated by larch, but it is comparable in k**kweight*chi*ftwin
        // The MSE is below 1.0e-4

        let larch_k_path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS_autobk_k_larch.txt";
        let larch_k = load_txt_f64(&larch_k_path, &PARAM_LOADTXT).unwrap();

        let k_expected = larch_k.get_col(0);
        let chi_expected = larch_k.get_col(1);

        let k = autobk.get_k().unwrap();
        let chi = autobk.get_chi_kweighted().unwrap();
        let ftwin = autobk.get_ftwin().unwrap();
        let kweight = autobk.get_kweight().unwrap();

        let chi_weighted = chi.component_mul(&ftwin);

        let chi_k2_weighted_expected = chi_expected
            .iter()
            .zip(k_expected.iter())
            .zip(ftwin.clone().iter())
            .map(|((x, y), z)| x * y.powi(*kweight) * z)
            .collect::<Vec<f64>>();

        let mse = chi_weighted
            .iter()
            .zip(chi_k2_weighted_expected.iter())
            .map(|(x, y)| (x - y).powi(2))
            .sum::<f64>()
            / chi_weighted.len() as f64;

        assert!(mse < CHI_MSE_TOL);
        Ok(())
    }

    #[test]
    fn test_autobk_knot_domain_bounds_and_ordering() {
        let spl_k = DVector::from_iterator(9, (0..9).map(|i| 1.5 + i as f64 * (11.0 / 8.0)));
        let order = 3usize;
        let knots = AUTOBK::build_knot_domain(&spl_k, order).unwrap();

        // Domain extension points are lower than qmin and upper than qmax.
        let qmin = spl_k[0];
        let qmax = spl_k[spl_k.len() - 1];
        assert!(knots[0] < qmin);
        assert!(knots[order - 1] < qmin);
        assert!(knots[order] >= qmin);
        assert!(*knots.last().unwrap() > qmax);

        // Knot vector must remain non-decreasing for spline construction.
        for i in 1..knots.len() {
            assert!(knots[i] >= knots[i - 1]);
        }

        // Interior knots stay inside [qmin, qmax].
        for knot in knots.iter().skip(order).take(spl_k.len() - order) {
            assert!(*knot >= qmin);
            assert!(*knot <= qmax);
        }
    }

    #[test]
    fn test_autobk_knot_domain_rejects_invalid_range() {
        let spl_k = DVector::from_vec(vec![2.0, 2.0, 2.0, 2.0, 2.0]);
        let err = AUTOBK::build_knot_domain(&spl_k, 3).unwrap_err();
        assert!(matches!(err, BackgroundError::SplineKnotsFailed { .. }));
    }

    #[test]
    fn test_autobk_automatic_knots_change_at_integer_boundary() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;
        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;
        let energy = spectrum.energy.as_ref().unwrap();
        let mu = spectrum.mu.as_ref().unwrap();

        // At rbkg=1, the ninth spline parameter is allowed only after
        // kmax crosses 4*pi. Compare actual fits to explicit knot counts,
        // including the common kmax=12 case that rounding handled differently.
        for (kmax, expected_knots) in [(12.0, 8), (12.565, 8), (12.568, 9)] {
            let mut automatic = AUTOBK::new();
            automatic.rbkg = Some(1.0);
            automatic.kmin = Some(0.0);
            automatic.kmax = Some(kmax);
            automatic.nclamp = Some(0);
            automatic.linear_fallback_to_lm = Some(false);
            let mut explicit = automatic.clone();
            explicit.nknots = Some(expected_knots);

            automatic.calc_background(energy, mu, &mut spectrum.normalization.clone())?;
            explicit.calc_background(energy, mu, &mut spectrum.normalization.clone())?;
            let actual = automatic.get_chi().unwrap();
            let expected = explicit.get_chi().unwrap();
            assert_eq!(actual.len(), expected.len());
            for (actual, expected) in actual.iter().zip(expected.iter()) {
                assert_abs_diff_eq!(actual, expected, epsilon = 1.0e-12);
            }
        }
        Ok(())
    }

    #[test]
    fn test_autobk_direct_solver_parity_with_legacy_lm() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;

        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();

        let mut legacy_norm = spectrum.normalization.clone();
        let mut direct_norm = spectrum.normalization.clone();

        let mut legacy = AUTOBK::new();
        legacy.solver = Some(AUTOBKSolver::LegacyLm);
        legacy.clamp_scale_policy = Some(AUTOBKClampScalePolicy::Fixed);
        legacy.calc_background(&energy, &mu, &mut legacy_norm)?;

        let mut direct = AUTOBK::new();
        direct.solver = Some(AUTOBKSolver::LinearDirect);
        direct.clamp_scale_policy = Some(AUTOBKClampScalePolicy::TwoPass);
        direct.calc_background(&energy, &mu, &mut direct_norm)?;

        let legacy_chi = legacy.get_chi().unwrap();
        let direct_chi = direct.get_chi().unwrap();
        assert_eq!(legacy_chi.len(), direct_chi.len());

        let mse = legacy_chi
            .iter()
            .zip(direct_chi.iter())
            .map(|(x, y)| (x - y).powi(2))
            .sum::<f64>()
            / legacy_chi.len() as f64;

        assert!(mse < CHI_MSE_TOL);
        Ok(())
    }

    #[cfg(feature = "trust-region")]
    #[test]
    fn test_autobk_trust_region_matches_legacy_lm() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;

        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();

        let mut legacy_norm = spectrum.normalization.clone();
        let mut trust_region_norm = spectrum.normalization.clone();

        let mut legacy = AUTOBK::new();
        legacy.solver = Some(AUTOBKSolver::LegacyLm);
        legacy.clamp_scale_policy = Some(AUTOBKClampScalePolicy::Fixed);
        legacy.calc_background(&energy, &mu, &mut legacy_norm)?;

        let mut trust_region = AUTOBK::new();
        trust_region.solver = Some(AUTOBKSolver::TrustRegionDogLeg);
        trust_region.clamp_scale_policy = Some(AUTOBKClampScalePolicy::Fixed);
        trust_region.calc_background(&energy, &mu, &mut trust_region_norm)?;

        let legacy_chi = legacy.get_chi().unwrap();
        let trust_region_chi = trust_region.get_chi().unwrap();
        assert_eq!(legacy_chi.len(), trust_region_chi.len());

        let mse = legacy_chi
            .iter()
            .zip(trust_region_chi.iter())
            .map(|(x, y)| (x - y).powi(2))
            .sum::<f64>()
            / legacy_chi.len() as f64;

        assert!(mse < CHI_MSE_TOL);
        Ok(())
    }

    #[cfg(feature = "trust-region")]
    #[test]
    fn test_autobk_direct_solver_fallback_to_trust_region() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;
        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();

        let mut normalization = spectrum.normalization.clone();

        let mut autobk = AUTOBK::new();
        autobk.clamp_scale_policy = Some(AUTOBKClampScalePolicy::Fixed);
        autobk.solver = Some(AUTOBKSolver::LinearDirect);
        autobk.linear_condition_limit = Some(1.0);
        autobk.linear_fallback_to_lm = Some(true);

        let result = autobk.calc_background(&energy, &mu, &mut normalization)?;
        assert!(result.get_chi().is_some());
        Ok(())
    }

    #[test]
    fn test_autobk_direct_solver_fallback_can_use_legacy_lm() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;
        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();
        let mut normalization = spectrum.normalization.clone();

        let mut autobk = AUTOBK::new();
        autobk.clamp_scale_policy = Some(AUTOBKClampScalePolicy::Fixed);
        autobk.solver = Some(AUTOBKSolver::LinearDirect);
        autobk.linear_condition_limit = Some(1.0);
        autobk.linear_fallback_to_lm = Some(true);
        autobk.linear_fallback_solver = Some(AUTOBKSolver::LegacyLm);

        let result = autobk.calc_background(&energy, &mu, &mut normalization)?;
        assert!(result.get_chi().is_some());
        Ok(())
    }

    #[test]
    fn test_autobk_deprecated_fallback_flag_preserves_legacy_lm_when_solver_missing() {
        let mut autobk = AUTOBK::new();
        autobk.linear_fallback_to_lm = Some(true);
        autobk.linear_fallback_solver = None;

        autobk.fill_parameter().unwrap();

        assert_eq!(autobk.linear_fallback_solver, Some(AUTOBKSolver::LegacyLm));
    }

    #[test]
    fn test_autobk_direct_solver_default_succeeds_without_fallback() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;
        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();
        let mut normalization = spectrum.normalization.clone();

        let mut autobk = AUTOBK::new();
        autobk.solver = Some(AUTOBKSolver::LinearDirect);
        autobk.linear_fallback_to_lm = Some(false);

        let result = autobk.calc_background(&energy, &mu, &mut normalization)?;
        assert!(result.get_chi().is_some());
        Ok(())
    }

    #[test]
    fn test_autobk_direct_solver_failure_without_fallback() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;
        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();
        let mut normalization = spectrum.normalization.clone();

        let mut autobk = AUTOBK::new();
        autobk.solver = Some(AUTOBKSolver::LinearDirect);
        autobk.linear_condition_limit = Some(1.0);
        autobk.linear_fallback_to_lm = Some(false);

        let err = autobk
            .calc_background(&energy, &mu, &mut normalization)
            .unwrap_err();
        assert!(matches!(
            err,
            BackgroundError::DirectSolverIllConditioned { .. }
                | BackgroundError::DirectSolverFailed { .. }
        ));
        Ok(())
    }

    #[test]
    fn test_autobk_direct_solver_deterministic_two_pass() -> Result<(), Box<dyn Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut spectrum = io::load_spectrum_QAS_trans(&path)?;
        spectrum
            .set_normalization_method(Some(normalization::NormalizationMethod::PrePostEdge(
                PrePostEdge::new(),
            )))?
            .normalize()?;

        let energy = spectrum.energy.clone().unwrap();
        let mu = spectrum.mu.clone().unwrap();

        let mut norm_a = spectrum.normalization.clone();
        let mut norm_b = spectrum.normalization.clone();

        let mut autobk_a = AUTOBK::new();
        autobk_a.solver = Some(AUTOBKSolver::LinearDirect);
        autobk_a.clamp_scale_policy = Some(AUTOBKClampScalePolicy::TwoPass);
        autobk_a.calc_background(&energy, &mu, &mut norm_a)?;

        let mut autobk_b = AUTOBK::new();
        autobk_b.solver = Some(AUTOBKSolver::LinearDirect);
        autobk_b.clamp_scale_policy = Some(AUTOBKClampScalePolicy::TwoPass);
        autobk_b.calc_background(&energy, &mu, &mut norm_b)?;

        let chi_a = autobk_a.get_chi().unwrap();
        let chi_b = autobk_b.get_chi().unwrap();
        assert_eq!(chi_a.len(), chi_b.len());

        chi_a.iter().zip(chi_b.iter()).for_each(|(a, b)| {
            assert_abs_diff_eq!(a, b, epsilon = 1.0e-12);
        });
        Ok(())
    }
}

#[cfg(test)]
#[path = "background/legacy_clamp_tests.rs"]
mod legacy_clamp_tests;