rexafs 0.2.5

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

use std::borrow::Borrow;
#[cfg_attr(debug_assertions, allow(dead_code, unused_imports))]
// Standard library dependencies
use std::error::Error;

// External dependencies
use easyfft::dyn_size::realfft::DynRealDft;
use nalgebra::DVector;
#[cfg(feature = "ndarray-compat")]
use ndarray::{ArrayBase, Ix1, ViewRepr};
use serde::{Deserialize, Serialize};

// load dependencies
use super::background;
use super::errors::{DataError, NormalizationError};
use super::io;
use super::lmutils;
use super::mathutils;
use super::normalization;
use super::nshare;
use super::tools;
use super::xafsutils;
use super::xrayfft;
use super::XAFSError;

// Load local traits
use mathutils::MathUtils;
use normalization::Normalization;

/// Data and processing parameters for a single XAS spectrum.
/// Also available as [`crate::Spectrum`]. Use [`Self::from_arrays`] for checked input.
/// The default is empty; processing creates missing default stage configurations.
/// Methods mutate this spectrum and getters never run a calculation implicitly.
///
/// Prefer getters over the legacy result fields: authoritative outputs live inside
/// `normalization`, `background`, `xftf` and `xftr`. Directly editing public inputs
/// or settings requires [`Self::invalidate_derived`] before processing again.
/// Raw arrays are a working baseline for interpolation, not an immutable archive:
/// calibration, deglitching, truncation, smoothing and rebinning can modify them.
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
#[serde(default)]
#[derive(Default)]
pub struct XASSpectrum {
    /// Optional display name; changing it does not affect numerical results.
    pub name: Option<String>,
    /// Owned baseline energy grid in eV. Data-treatment methods may modify it.
    pub raw_energy: Option<DVector<f64>>,
    /// Owned baseline absorption values, paired with `raw_energy`.
    pub raw_mu: Option<DVector<f64>>,
    /// Owned current energy grid in eV used by processing stages.
    pub energy: Option<DVector<f64>>,
    /// Owned current absorption values, paired with `energy`; units match the input.
    pub mu: Option<DVector<f64>>,
    /// Selected or estimated edge energy in eV. Prefer [`Self::set_e0`] to edit it.
    pub e0: Option<f64>,
    /// Legacy result slot; use [`Self::k`] for the background wave-number grid in Å⁻¹.
    pub k: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::chi`] for dimensionless unweighted EXAFS.
    pub chi: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::chi_kweighted`] to calculate weighted EXAFS.
    pub chi_kweighted: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::chir`] for the stored complex Fourier data.
    pub chi_r: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::chir_mag`] for Fourier magnitudes.
    pub chi_r_mag: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::chir_real`] for the real Fourier component.
    pub chi_r_re: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::chir_imag`] for the imaginary Fourier component.
    pub chi_r_im: Option<DVector<f64>>,
    /// Legacy result slot; use [`Self::q`] for the inverse-transform grid in Å⁻¹.
    pub q: Option<DVector<f64>>,
    /// Normalization settings and cached outputs; `None` selects default pre/post-edge normalization when needed.
    pub normalization: Option<normalization::NormalizationMethod>,
    /// Background settings and cached outputs; `None` selects default AUTOBK when needed.
    pub background: Option<background::BackgroundMethod>,
    /// Forward-transform settings and outputs; `None` selects [`xrayfft::XrayFFTF::default`].
    pub xftf: Option<xrayfft::XrayFFTF>,
    /// Inverse-transform settings and outputs; `None` selects [`xrayfft::XrayFFTR::default`].
    pub xftr: Option<xrayfft::XrayFFTR>,
    /// Accumulated energy shift (eV) applied by `shift_energy`/`calibrate`/`align_to`.
    pub energy_shift: f64,
    /// Per-point spread stored by merge/rebin, in input absorption units.
    /// Other data edits do not consistently propagate, resize or clear this field;
    /// verify alignment with `energy`/`mu` before reuse. It is not automatically
    /// consumed as an uncertainty model by spectrum processing.
    pub mu_stddev: Option<DVector<f64>>,
    /// Marker set by rebinning. Later data replacement does not reset it, so it
    /// records a past operation rather than validating the current arrays.
    pub rebinned: bool,
    /// Explicit normalization scale, distinct from the scale inferred by a stage.
    #[serde(skip_serializing_if = "Option::is_none")]
    normalization_edge_step_override: Option<f64>,
    /// Last normalization output, used to recognize later public-field edits.
    #[serde(skip_serializing_if = "Option::is_none")]
    normalization_edge_step_last_result: Option<f64>,
}

impl XASSpectrum {
    fn validate_energy_mu_inputs(
        energy: &DVector<f64>,
        mu: &DVector<f64>,
    ) -> Result<(), XAFSError> {
        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(())
    }

    /// Create an empty spectrum with no data, selected methods or calculated results.
    /// Use [`Self::from_arrays`] to construct checked input in one call.
    pub fn new() -> XASSpectrum {
        XASSpectrum::default()
    }

    /// Create a spectrum from finite, equal-length arrays with strictly increasing
    /// energy in eV. Unlike the legacy setter, this validates before indexing.
    pub fn from_arrays(energy: &[f64], mu: &[f64]) -> Result<Self, XAFSError> {
        let energy = DVector::from_column_slice(energy);
        let mu = DVector::from_column_slice(mu);
        Self::validate_energy_mu_inputs(&energy, &mu)?;
        for index in 1..energy.len() {
            if energy[index] == energy[index - 1] {
                return Err(DataError::DuplicateEnergy {
                    index,
                    energy: energy[index],
                }
                .into());
            }
        }
        let mut spectrum = Self::new();
        spectrum.set_spectrum(energy, mu);
        Ok(spectrum)
    }

    /// Borrow the background k grid without cloning its buffer (Å⁻¹).
    /// Returns `None` before a valid AUTOBK result; no calculation is triggered.
    pub fn k(&self) -> Option<&[f64]> {
        let background::BackgroundMethod::AUTOBK(autobk) = self.background.as_ref()? else {
            return None;
        };
        #[cfg(not(feature = "ndarray-compat"))]
        {
            Some(autobk.k.as_ref()?.as_slice())
        }
        #[cfg(feature = "ndarray-compat")]
        {
            autobk.k.as_ref()?.as_slice()
        }
    }

    /// Borrow dimensionless, unweighted χ(k) without cloning its buffer.
    /// Returns `None` before a valid AUTOBK result; no calculation is triggered.
    pub fn chi(&self) -> Option<&[f64]> {
        let background::BackgroundMethod::AUTOBK(autobk) = self.background.as_ref()? else {
            return None;
        };
        #[cfg(not(feature = "ndarray-compat"))]
        {
            Some(autobk.chi.as_ref()?.as_slice())
        }
        #[cfg(feature = "ndarray-compat")]
        {
            autobk.chi.as_ref()?.as_slice()
        }
    }

    /// Set the display name without invalidating numerical results.
    pub fn set_name<S: Into<String>>(&mut self, name: S) -> &mut Self {
        self.name = Some(name.into());
        self
    }

    /// Replace the baseline and working arrays, sorting energy and absorption together.
    /// Energy is in eV. This legacy setter takes ownership after conversion and clones
    /// the baseline into working arrays. It does not check lengths or finite values;
    /// prefer [`Self::from_arrays`] for checked input. Clears E0 and derived results
    /// while retaining other stage settings.
    ///
    /// # Panics
    /// May panic when sorting mismatched arrays. Supply paired finite arrays.
    pub fn set_spectrum<T: Into<DVector<f64>>, M: Into<DVector<f64>>>(
        &mut self,
        energy: T,
        mu: M,
    ) -> &mut Self {
        let raw_energy = energy.into();
        let raw_mu = mu.into();

        if !raw_energy.is_sorted() {
            let sort_idx = raw_energy.argsort();
            // For DVector, we need to manually sort by indices
            self.raw_energy = Some(DVector::from_iterator(
                sort_idx.len(),
                sort_idx.iter().map(|&i| raw_energy[i]),
            ));
            self.raw_mu = Some(DVector::from_iterator(
                sort_idx.len(),
                sort_idx.iter().map(|&i| raw_mu[i]),
            ));
        } else {
            self.raw_energy = Some(raw_energy);
            self.raw_mu = Some(raw_mu);
        }
        self.energy = self.raw_energy.clone();
        self.mu = self.raw_mu.clone();
        self.e0 = None;
        if let Some(method) = self.normalization.as_mut() {
            method.set_e0(None);
        }
        if let Some(background::BackgroundMethod::AUTOBK(method)) = self.background.as_mut() {
            method.ek0 = None;
        }
        self.invalidate_derived()
    }

    /// Linearly interpolate baseline absorption onto an owned energy grid in eV.
    /// The baseline arrays remain unchanged; successful interpolation invalidates
    /// derived results. The interpolation helper holds endpoint values outside the
    /// baseline range. Missing baseline data or interpolation failures return an error.
    /// The working energy grid is assigned before interpolation, even if it fails.
    pub fn interpolate_spectrum<T: Into<DVector<f64>>>(
        &mut self,
        energy: T,
    ) -> Result<&mut Self, XAFSError> {
        self.energy = Some(energy.into());
        let energy = self.energy.as_ref().ok_or_else(|| DataError::MissingData {
            field: "energy".to_string(),
        })?;
        let mu = self.raw_mu.as_ref().ok_or_else(|| DataError::MissingData {
            field: "raw_mu".to_string(),
        })?;
        let knot = self
            .raw_energy
            .as_ref()
            .ok_or_else(|| DataError::MissingData {
                field: "raw_energy".to_string(),
            })?;

        let interpolated = energy
            .interpolate(knot.as_slice(), mu.as_slice())
            .map_err(|e| super::errors::MathError::SplineEvalFailed {
                x: 0.0,
                reason: e.to_string(),
            })?;
        self.mu = Some(interpolated);

        Ok(self.invalidate_derived())
    }

    /// Set the edge energy in eV and propagate it into existing stage configurations.
    /// Invalidates normalization, background and Fourier results. The value is stored
    /// as supplied; normalization subsequently rejects a non-finite E0 or one outside
    /// the measured range. Use [`Self::find_e0`] for an automatic estimate.
    pub fn set_e0<S: Into<f64>>(&mut self, e0: S) -> &mut Self {
        self.invalidate_derived();
        self.e0 = Some(e0.into());
        if let Some(method) = self.normalization.as_mut() {
            method.set_e0(self.e0);
        }
        if let Some(background::BackgroundMethod::AUTOBK(method)) = self.background.as_mut() {
            method.ek0 = self.e0;
        }
        self
    }

    /// Estimate E0 in eV from the current spectrum using the core edge detector.
    /// Propagates the result through [`Self::set_e0`] and invalidates derived results.
    /// Missing, non-finite, mismatched or decreasing data return a typed error; a usable
    /// absorption edge and enough points are also required by the detector.
    pub fn find_e0(&mut self) -> Result<&mut Self, XAFSError> {
        let energy = self.energy.as_ref().ok_or_else(|| DataError::MissingData {
            field: "energy".to_string(),
        })?;
        let mu = self.mu.as_ref().ok_or_else(|| DataError::MissingData {
            field: "mu".to_string(),
        })?;
        Self::validate_energy_mu_inputs(energy, mu)?;
        let e0 = xafsutils::find_e0(energy, mu)?;
        Ok(self.set_e0(e0))
    }

    fn find_energy_step(
        &mut self,
        frac_ignore: Option<f64>,
        nave: Option<usize>,
    ) -> Result<f64, XAFSError> {
        let energy = self.energy.as_ref().ok_or_else(|| DataError::MissingData {
            field: "energy".to_string(),
        })?;
        if energy.len() < 2 {
            return Err(DataError::InsufficientData {
                min: 2,
                actual: energy.len(),
            }
            .into());
        }
        let non_finite = energy
            .iter()
            .enumerate()
            .filter_map(|(index, value)| (!value.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(xafsutils::find_energy_step(energy, frac_ignore, nave, None))
    }

    /// Take ownership of normalization settings and clear dependent results.
    /// Pass [`crate::PrePostEdge`] directly, a method enum, or an optional enum.
    /// `None` selects default pre/post-edge normalization. A configured edge energy
    /// takes precedence over the existing E0; an explicit edge step remains an override.
    /// No normalization is performed by this setter.
    pub fn set_normalization_method(
        &mut self,
        method: impl Into<Option<normalization::NormalizationMethod>>,
    ) -> Result<&mut Self, XAFSError> {
        let method = method.into();
        self.invalidate_derived();
        if let Some(method) = method {
            self.normalization = Some(method);
        } else {
            let normalization_method = normalization::PrePostEdge::new();
            self.normalization = Some(normalization::NormalizationMethod::PrePostEdge(
                normalization_method,
            ));
        }

        self.normalization_edge_step_override =
            self.normalization.as_ref().and_then(|m| m.get_edge_step());
        let e0 = self
            .normalization
            .as_ref()
            .and_then(|m| m.get_e0())
            .or(self.e0);
        if e0 != self.e0 {
            if let Some(background::BackgroundMethod::AUTOBK(method)) = self.background.as_mut() {
                method.ek0 = e0;
            }
        }
        self.e0 = e0;
        if let Some(normalization_method) = self.normalization.as_mut() {
            if e0.is_some() {
                normalization_method.set_e0(e0);
            }
        } else {
            return Err(DataError::MissingData {
                field: "normalization method".to_string(),
            }
            .into());
        }

        Ok(self)
    }

    /// Normalize using the selected method, resolving E0 and defaults as needed.
    /// The default fits pre/post-edge curves and expresses absorption in edge-step
    /// units. Call [`Self::norm`] or [`Self::flat`] for owned outputs after success.
    /// Recomputing invalidates background and Fourier results. Missing/invalid data,
    /// an unusable fitting range, or an unsupported method returns a typed error.
    pub fn normalize(&mut self) -> Result<&mut Self, XAFSError> {
        // Capture explicitly configured edge_step before the algorithm fills it.
        if let Some(method) = &self.normalization {
            if method.get_norm().is_none() {
                self.normalization_edge_step_override = method.get_edge_step();
            }
        }
        self.invalidate_derived();
        if self.normalization.is_none() {
            self.set_normalization_method(None)?;
        }

        let configured_e0 = self.normalization.as_ref().and_then(|m| m.get_e0());
        // Let the selected normalization method estimate an unset E0. In
        // ndarray compatibility mode this preserves its legacy edge detector.

        let energy = self.energy.as_ref().ok_or_else(|| DataError::MissingData {
            field: "energy".to_string(),
        })?;
        let mu = self.mu.as_ref().ok_or_else(|| DataError::MissingData {
            field: "mu".to_string(),
        })?;
        Self::validate_energy_mu_inputs(energy, mu)?;

        if let Some(e0) = self.e0.or(configured_e0) {
            let data_min = energy[0];
            let data_max = energy[energy.len() - 1];
            if !e0.is_finite() || e0 <= data_min || e0 >= data_max {
                return Err(NormalizationError::E0OutOfRange {
                    e0,
                    data_min,
                    data_max,
                }
                .into());
            }
        }

        let mut method = self
            .normalization
            .clone()
            .ok_or_else(|| DataError::MissingData {
                field: "normalization method".to_string(),
            })?;
        if let Err(error) = method.normalize(energy, mu) {
            self.invalidate_derived();
            return Err(error.into());
        }
        self.normalization = Some(method);
        self.e0 = self.normalization.as_ref().and_then(|m| m.get_e0());
        self.normalization_edge_step_last_result =
            self.normalization.as_ref().and_then(|m| m.get_edge_step());
        Ok(self)
    }

    /// Take ownership of background settings and clear background and Fourier results.
    /// Pass [`crate::AUTOBK`] directly, a method enum, or an optional enum.
    /// `None` selects default AUTOBK; existing normalization is retained. No background
    /// calculation is performed by this setter.
    pub fn set_background_method(
        &mut self,
        method: impl Into<Option<background::BackgroundMethod>>,
    ) -> Result<&mut Self, XAFSError> {
        let method = method.into();
        self.invalidate_background();
        if let Some(method) = method {
            self.background = Some(method);
        } else {
            let backgound_method = background::AUTOBK::new();
            self.background = Some(background::BackgroundMethod::AUTOBK(backgound_method));
        }

        Ok(self)
    }

    /// Remove the smooth background using the selected method; normalize if needed.
    /// Default nalgebra AUTOBK minimizes low-R content with a fixed endpoint penalty;
    /// the optional `ndarray-compat` backend retains its historical clamp model.
    /// Successful results are dimensionless unweighted [`Self::chi`] on [`Self::k`].
    /// Recomputes this stage and clears Fourier results on every call. Missing data,
    /// insufficient coverage, invalid settings or solver failure return a typed error.
    pub fn calc_background(&mut self) -> Result<&mut Self, XAFSError> {
        self.invalidate_background();
        if self
            .normalization
            .as_ref()
            .and_then(|m| m.get_norm())
            .is_none()
        {
            self.normalize()?;
        }
        if self.background.is_none() {
            self.set_background_method(None)?;
        }

        let energy = self.energy.as_ref().ok_or_else(|| DataError::MissingData {
            field: "energy".to_string(),
        })?;
        let mu = self.mu.as_ref().ok_or_else(|| DataError::MissingData {
            field: "mu".to_string(),
        })?;
        Self::validate_energy_mu_inputs(energy, mu)?;

        self.background
            .as_mut()
            .ok_or_else(|| DataError::MissingData {
                field: "background method".to_string(),
            })?
            .calc_background(energy, mu, &mut self.normalization)?;

        Ok(self)
    }

    /// Take ownership of forward settings and clear forward/inverse results.
    /// Existing normalization and background results are retained. Values already
    /// inferred inside `parameters`, such as `kstep`, stay explicit; use a fresh
    /// configuration or reset such fields to `None` to request new inference.
    pub fn set_fft(&mut self, parameters: xrayfft::XrayFFTF) -> &mut Self {
        self.xftf = Some(parameters);
        self.invalidate_fft();
        self
    }

    /// Weight and window χ(k), then transform it to complex Fourier distance R.
    /// Computes missing normalization/background stages first, using their defaults.
    /// The unnormalized negative-exponent FFT is multiplied by `kstep / sqrt(pi)`;
    /// no extra division by the FFT length or window area is applied.
    /// Default settings use k-weight 2, a Kaiser–Bessel window and 2048 samples.
    /// Invalid input/settings or a failed prerequisite returns a typed error.
    /// Every call recomputes the forward transform and clears inverse results.
    /// Successful calls retain resolved settings: an inferred `kstep` is reused
    /// after later background-grid changes unless reset through [`Self::set_fft`].
    pub fn fft(&mut self) -> Result<&mut Self, XAFSError> {
        self.invalidate_fft();
        if self.k().is_none() || self.chi().is_none() {
            self.calc_background()?;
        }
        // Work on a copy so errors retain the caller's transform parameters.
        let mut xftf = self.xftf.clone().unwrap_or_default();

        #[cfg(feature = "ndarray-compat")]
        {
            let k = self.k_view().ok_or_else(|| DataError::MissingData {
                field: "k (need to calculate background first)".to_string(),
            })?;
            let chi = self.chi_view().ok_or_else(|| DataError::MissingData {
                field: "chi (need to calculate background first)".to_string(),
            })?;
            xftf.xftf(k, chi)?;
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            let k = self.k().ok_or_else(|| DataError::MissingData {
                field: "k (need to calculate background first)".to_string(),
            })?;
            let chi = self.chi().ok_or_else(|| DataError::MissingData {
                field: "chi (need to calculate background first)".to_string(),
            })?;
            xftf.xftf(
                &DVector::from_column_slice(k),
                &DVector::from_column_slice(chi),
            )?;
        }

        self.xftf = Some(xftf);
        Ok(self)
    }

    /// Take ownership of inverse settings, preserving forward results and clearing q/chi(q).
    /// Previously inferred `kstep`/`nfft` values in the supplied settings remain
    /// explicit. Pass a fresh configuration after changing forward-grid geometry
    /// when you want the inverse grid to be inferred again.
    pub fn set_ifft(&mut self, mut parameters: xrayfft::XrayFFTR) -> &mut Self {
        parameters.q = None;
        parameters.chiq = None;
        parameters.rwin = None;
        self.q = None;
        self.xftr = Some(parameters);
        self
    }

    /// Filter χ(R) and return a real inverse transform, computing missing stages first.
    /// Uses [`xrayfft::XrayFFTR`] defaults unless inverse settings are configured.
    /// This preserves the forward weighting/window and is not an unweighted χ(k)
    /// reconstruction. Invalid inverse settings or prerequisites return an error.
    /// Resolved inverse grid settings persist across calls; after changing forward
    /// `nfft` or spacing, reset them with [`Self::set_ifft`] to request inference.
    pub fn ifft(&mut self) -> Result<&mut Self, XAFSError> {
        if self.chir().is_none() {
            self.fft()?;
        }

        let xftf = self.xftf.as_ref().ok_or_else(|| DataError::MissingData {
            field: "xftf configuration".to_string(),
        })?;

        if self.xftr.is_none() {
            self.xftr = Some(xrayfft::XrayFFTR::new());
        }

        #[cfg(feature = "ndarray-compat")]
        self.xftr
            .as_mut()
            .ok_or_else(|| DataError::MissingData {
                field: "xftr configuration".to_string(),
            })?
            .xftr(
                xftf.get_r().ok_or_else(|| DataError::MissingData {
                    field: "r (fft() may have failed)".to_string(),
                })?,
                xftf.get_chir().ok_or_else(|| DataError::MissingData {
                    field: "chi_r (fft() may have failed)".to_string(),
                })?,
            )?;

        #[cfg(not(feature = "ndarray-compat"))]
        self.xftr
            .as_mut()
            .ok_or_else(|| DataError::MissingData {
                field: "xftr configuration".to_string(),
            })?
            .xftr(
                xftf.get_r().ok_or_else(|| DataError::MissingData {
                    field: "r (fft() may have failed)".to_string(),
                })?,
                xftf.get_chir().ok_or_else(|| DataError::MissingData {
                    field: "chi_r (fft() may have failed)".to_string(),
                })?,
            )?;

        Ok(self)
    }

    // -----------------------------------------------------------------------
    // Athena-style data-processing tools (see `xafs::tools`)
    // -----------------------------------------------------------------------

    /// Clear every result derived from `energy`/`mu` (normalization outputs,
    /// background, χ(k), χ(R), χ(q)) while keeping the stage parameters, so
    /// the pipeline recomputes from the modified data. Resolved automatic settings
    /// are retained along with explicit settings: clearing a result does not set
    /// its inferred `kstep`, `nfft` or fitting ranges back to `None`. This also leaves
    /// `mu_stddev`, `rebinned` and the selected E0 unchanged.
    pub fn invalidate_derived(&mut self) -> &mut Self {
        if let Some(last_result) = self.normalization_edge_step_last_result.take() {
            let current = self.normalization.as_ref().and_then(|m| m.get_edge_step());
            if current != Some(last_result) {
                // The caller changed (or cleared) the public parameter after
                // normalization. Preserve that choice instead of the old scale.
                self.normalization_edge_step_override = current;
            }
        }
        match self.normalization.as_mut() {
            Some(normalization::NormalizationMethod::PrePostEdge(p)) => {
                if p.norm.is_some() {
                    p.edge_step = self.normalization_edge_step_override;
                }
                p.pre_edge = None;
                p.post_edge = None;
                p.norm = None;
                p.flat = None;
                p.pre_coefficients = None;
                p.norm_coefficients = None;
            }
            Some(normalization::NormalizationMethod::MBack(m)) => {
                if m.norm.is_some() {
                    m.edge_step = self.normalization_edge_step_override;
                }
                m.norm = None;
                m.flat = None;
            }
            None => {}
        }
        self.invalidate_background();
        self
    }

    fn invalidate_background(&mut self) {
        self.k = None;
        self.chi = None;
        self.chi_kweighted = None;
        if let Some(background::BackgroundMethod::AUTOBK(a)) = self.background.as_mut() {
            a.bkg = None;
            a.chie = None;
            a.k = None;
            a.chi = None;
        }
        self.invalidate_fft();
    }

    fn invalidate_fft(&mut self) {
        self.chi_r = None;
        self.chi_r_mag = None;
        self.chi_r_re = None;
        self.chi_r_im = None;
        self.q = None;
        if let Some(f) = self.xftf.as_mut() {
            f.r = None;
            f.chir = None;
            f.chir_mag = None;
            f.kwin = None;
        }
        if let Some(r) = self.xftr.as_mut() {
            r.q = None;
            r.chiq = None;
            r.rwin = None;
        }
    }

    fn working_pair(&self) -> Result<(&DVector<f64>, &DVector<f64>), XAFSError> {
        let energy = self.energy.as_ref().ok_or_else(|| DataError::MissingData {
            field: "energy".to_string(),
        })?;
        let mu = self.mu.as_ref().ok_or_else(|| DataError::MissingData {
            field: "mu".to_string(),
        })?;
        Self::validate_energy_mu_inputs(energy, mu)?;
        Ok((energy, mu))
    }

    fn raw_grid_matches_working(&self) -> bool {
        match (&self.raw_energy, &self.energy) {
            (Some(raw), Some(energy)) => raw == energy,
            _ => false,
        }
    }

    /// Shift the energy axis (working and raw) by `delta_ev`, moving `e0` and
    /// the e0-like stage parameters along with it. The shift accumulates in
    /// `energy_shift`. Derived results are invalidated. The shift is in eV and
    /// must be finite; this setter stores it without immediate validation.
    pub fn shift_energy(&mut self, delta_ev: f64) -> &mut Self {
        if let Some(e) = self.energy.as_mut() {
            e.add_scalar_mut(delta_ev);
        }
        if let Some(e) = self.raw_energy.as_mut() {
            e.add_scalar_mut(delta_ev);
        }
        self.energy_shift += delta_ev;
        if let Some(e0) = self.e0.as_mut() {
            *e0 += delta_ev;
        }
        if let Some(norm) = self.normalization.as_mut() {
            let e0 = norm.get_e0().map(|e0| e0 + delta_ev);
            norm.set_e0(e0);
        }
        if let Some(background::BackgroundMethod::AUTOBK(a)) = self.background.as_mut() {
            a.ek0 = a.ek0.map(|e| e + delta_ev);
        }
        self.invalidate_derived()
    }

    /// Energy of the requested edge feature on the current `energy`/`mu`.
    /// `HalfStep` normalizes a copy of the spectrum if no flattened μ is available.
    pub fn edge_feature_energy(&self, feature: tools::EdgeFeature) -> Result<f64, XAFSError> {
        let (energy, mu) = self.working_pair()?;
        match feature {
            tools::EdgeFeature::DerivativeMax => tools::derivative_max_energy(energy, mu),
            tools::EdgeFeature::SecondDerivativeZero => {
                tools::second_derivative_zero_energy(energy, mu)
            }
            tools::EdgeFeature::HalfStep => {
                let e0 = match self.e0 {
                    Some(e0) => e0,
                    None => tools::derivative_max_energy(energy, mu)?,
                };
                let flat = match self.flat() {
                    Some(flat) if flat.len() == energy.len() => flat,
                    _ => {
                        let mut tmp = self.clone();
                        tmp.invalidate_derived();
                        tmp.normalize()?;
                        tmp.flat().ok_or_else(|| DataError::MissingData {
                            field: "flat".to_string(),
                        })?
                    }
                };
                tools::half_step_energy(energy, &flat, e0)
            }
        }
    }

    /// Calibrate: shift the spectrum so that `feature` lands on `target_ev`
    /// and set the spectrum/normalization E0 to `target_ev`, in eV. Returns the
    /// shift applied in eV and invalidates derived results. A feature-detection
    /// failure is returned before shifting. The target must be finite; no immediate
    /// target validation is performed. Existing AUTOBK `ek0` is shifted with the
    /// axes and can differ from the target when its old reference differed from
    /// the selected feature; use [`Self::set_e0`] to synchronize it explicitly.
    pub fn calibrate(
        &mut self,
        feature: tools::EdgeFeature,
        target_ev: f64,
    ) -> Result<f64, XAFSError> {
        let current = self.edge_feature_energy(feature)?;
        let shift = target_ev - current;
        self.shift_energy(shift);
        self.e0 = Some(target_ev);
        if let Some(norm) = self.normalization.as_mut() {
            norm.set_e0(Some(target_ev));
        }
        Ok(shift)
    }

    /// Align this spectrum to `reference` by overlaying dμ/dE within
    /// `window` (eV, relative to the reference e0). The best shift (searched
    /// on a coarse ±20 eV interval with 0.1 eV steps) is applied with
    /// [`Self::shift_energy`] and returned in eV. Refinement can move slightly
    /// outside that interval. The reference is not mutated; insufficient coverage
    /// or an invalid prerequisite returns an error. See [`tools::find_energy_shift`]
    /// for the sign convention and free amplitude scaling.
    pub fn align_to(
        &mut self,
        reference: &XASSpectrum,
        window: (f64, f64),
    ) -> Result<f64, XAFSError> {
        let (e_dat, mu_dat) = self.working_pair()?;
        let (e_ref, mu_ref) = reference.working_pair()?;
        let ref_e0 = match reference.e0 {
            Some(e0) => e0,
            None => tools::derivative_max_energy(e_ref, mu_ref)?,
        };
        let shift = tools::find_energy_shift(
            e_dat,
            &tools::dmude(e_dat, mu_dat),
            e_ref,
            &tools::dmude(e_ref, mu_ref),
            ref_e0 + window.0,
            ref_e0 + window.1,
            20.0,
            0.1,
        )?;
        self.shift_energy(shift);
        Ok(shift)
    }

    /// Remove the points nearest to `energies` from the working and raw
    /// arrays. Returns the number of working points removed.
    fn remove_points_at(&mut self, energies: &[f64]) -> Result<usize, XAFSError> {
        if energies.is_empty() {
            return Ok(0);
        }
        let (energy, mu) = self.working_pair()?;
        let idx = tools::nearest_indices(energy, energies);
        if energy.len() - idx.len() < 2 {
            return Err(DataError::InsufficientData {
                min: 2,
                actual: energy.len() - idx.len(),
            }
            .into());
        }
        let new_energy = tools::remove_indices(energy, &idx);
        let new_mu = tools::remove_indices(mu, &idx);
        let new_raw = match (self.raw_energy.as_ref(), self.raw_mu.as_ref()) {
            (Some(raw_e), Some(raw_mu))
                if raw_e.len() == raw_mu.len() && raw_e.len() > idx.len() + 1 =>
            {
                let raw_idx = tools::nearest_indices(raw_e, energies);
                Some((
                    tools::remove_indices(raw_e, &raw_idx),
                    tools::remove_indices(raw_mu, &raw_idx),
                ))
            }
            _ => None,
        };
        if let Some((re, rm)) = new_raw {
            self.raw_energy = Some(re);
            self.raw_mu = Some(rm);
        }
        self.energy = Some(new_energy);
        self.mu = Some(new_mu);
        self.invalidate_derived();
        Ok(idx.len())
    }

    /// Remove working samples nearest to target energies in eV.
    /// Returns the number of distinct working samples removed; repeated targets
    /// remove a sample once and out-of-range targets select an endpoint. When
    /// usable raw arrays exist, removes their nearest samples too. Invalidates
    /// derived results and errors if fewer than two working samples would remain.
    /// Stored `mu_stddev` is not resized or recalculated.
    pub fn deglitch_points(&mut self, energies_to_remove: &[f64]) -> Result<usize, XAFSError> {
        self.remove_points_at(energies_to_remove)
    }

    /// Remove every working point in an inclusive energy interval, in eV.
    /// Reversed bounds are swapped. Delegates the selected energies to
    /// [`Self::deglitch_points`], including its raw-array and uncertainty behavior.
    pub fn deglitch_range(&mut self, e_lo: f64, e_hi: f64) -> Result<usize, XAFSError> {
        let (energy, _) = self.working_pair()?;
        let targets: Vec<f64> = tools::indices_in_range(energy, e_lo, e_hi)
            .into_iter()
            .map(|i| energy[i])
            .collect();
        self.remove_points_at(&targets)
    }

    /// Athena's margin deglitch: fit a line to μ(E) over `[e_lo, e_hi]` and
    /// remove the points lying more than `upper_margin` above or
    /// `lower_margin` below it. Bounds are in eV and margins in absorption units;
    /// negative margins use their absolute values. Returns removed energies in eV.
    /// At least two selected points must support the fitted line and two working
    /// samples must remain. Data mutation follows [`Self::deglitch_points`].
    pub fn deglitch_margin(
        &mut self,
        e_lo: f64,
        e_hi: f64,
        upper_margin: f64,
        lower_margin: f64,
    ) -> Result<Vec<f64>, XAFSError> {
        let (energy, mu) = self.working_pair()?;
        let removed: Vec<f64> =
            tools::margin_outliers(energy, mu, e_lo, e_hi, upper_margin, lower_margin)?
                .into_iter()
                .map(|i| energy[i])
                .collect();
        self.remove_points_at(&removed)?;
        Ok(removed)
    }

    /// Truncate: keep only points with `before <= E <= after` (either bound
    /// may be `None` for no bound). Bounds are in eV and are not automatically
    /// swapped. Errors before mutation if fewer than two working samples remain.
    /// Raw arrays are truncated too when at least two raw points remain; otherwise
    /// they stay unchanged. Clears derived results but does not resize `mu_stddev`.
    pub fn truncate(
        &mut self,
        before: Option<f64>,
        after: Option<f64>,
    ) -> Result<&mut Self, XAFSError> {
        let lo = before.unwrap_or(f64::NEG_INFINITY);
        let hi = after.unwrap_or(f64::INFINITY);
        let keep = |e: &DVector<f64>, m: &DVector<f64>| -> (DVector<f64>, DVector<f64>) {
            let idx: Vec<usize> = e
                .iter()
                .enumerate()
                .filter(|(_, v)| **v >= lo && **v <= hi)
                .map(|(i, _)| i)
                .collect();
            (
                DVector::from_iterator(idx.len(), idx.iter().map(|&i| e[i])),
                DVector::from_iterator(idx.len(), idx.iter().map(|&i| m[i])),
            )
        };
        let (energy, mu) = self.working_pair()?;
        let (new_energy, new_mu) = keep(energy, mu);
        if new_energy.len() < 2 {
            return Err(DataError::InsufficientData {
                min: 2,
                actual: new_energy.len(),
            }
            .into());
        }
        let new_raw = match (self.raw_energy.as_ref(), self.raw_mu.as_ref()) {
            (Some(raw_e), Some(raw_mu)) if raw_e.len() == raw_mu.len() => Some(keep(raw_e, raw_mu)),
            _ => None,
        };
        if let Some((re, rm)) = new_raw {
            if re.len() >= 2 {
                self.raw_energy = Some(re);
                self.raw_mu = Some(rm);
            }
        }
        self.energy = Some(new_energy);
        self.mu = Some(new_mu);
        Ok(self.invalidate_derived())
    }

    /// Rebin onto Athena's three-region grid (see [`tools::rebin`]). The raw
    /// arrays are replaced by the rebinned data, `mu_stddev` holds the
    /// per-bin standard deviation and `rebinned` is set. E0 comes from `cfg.e0`,
    /// then the current spectrum E0, then an automatic estimate used for the grid.
    /// The automatically estimated grid E0 is not copied into `self.e0` by this
    /// method. Errors from [`tools::rebin`] leave the arrays unchanged. Subsequent
    /// stages are invalidated, while their resolved parameters remain stored.
    pub fn rebin(&mut self, cfg: &tools::RebinConfig) -> Result<&mut Self, XAFSError> {
        let (energy, mu) = self.working_pair()?;
        let cfg = tools::RebinConfig {
            e0: cfg.e0.or(self.e0),
            ..*cfg
        };
        let out = tools::rebin(energy, mu, &cfg)?;
        self.raw_energy = Some(out.energy.clone());
        self.raw_mu = Some(out.mu.clone());
        self.energy = Some(out.energy);
        self.mu = Some(out.mu);
        self.mu_stddev = Some(out.stddev);
        self.rebinned = true;
        Ok(self.invalidate_derived())
    }

    /// Clone the spectrum and apply [`Self::rebin`] to the clone.
    /// A named result receives the suffix ` (rebinned)`; the original arrays,
    /// settings and results remain unchanged. Returns the same rebin errors.
    pub fn rebinned(&self, cfg: &tools::RebinConfig) -> Result<XASSpectrum, XAFSError> {
        let mut out = self.clone();
        out.rebin(cfg)?;
        if let Some(name) = self.name.as_deref() {
            out.set_name(format!("{name} (rebinned)"));
        }
        Ok(out)
    }

    /// Replace working μ(E) with a Lorentzian/Gaussian/Voigt smoothed copy.
    /// Width definitions and defaults (`sigma = 1` eV, `gamma = sigma`) follow
    /// [`tools::smooth_mu`]. Matching raw arrays receive the same result; a
    /// different usable raw grid is smoothed independently, otherwise raw data is
    /// retained. Invalidates derived results. A smoothing failure returns before
    /// arrays are replaced. Stored `mu_stddev` is not propagated through the filter.
    pub fn smooth_mu(
        &mut self,
        form: xafsutils::ConvolveForm,
        sigma: Option<f64>,
        gamma: Option<f64>,
    ) -> Result<&mut Self, XAFSError> {
        let (energy, mu) = self.working_pair()?;
        let smoothed = tools::smooth_mu(energy, mu, form, sigma, gamma)?;
        let new_raw_mu = if self.raw_grid_matches_working() {
            Some(smoothed.clone())
        } else {
            match (self.raw_energy.as_ref(), self.raw_mu.as_ref()) {
                (Some(raw_e), Some(raw_mu)) if raw_e.len() == raw_mu.len() && raw_e.len() >= 3 => {
                    Some(tools::smooth_mu(raw_e, raw_mu, form, sigma, gamma)?)
                }
                _ => None,
            }
        };
        if new_raw_mu.is_some() {
            self.raw_mu = new_raw_mu;
        }
        self.mu = Some(smoothed);
        Ok(self.invalidate_derived())
    }

    /// Return the selected or estimated edge energy in eV, or `None` before resolution.
    pub fn e0(&self) -> Option<f64> {
        self.e0
    }

    /// Copy normalized absorption in edge-step units on the current energy grid.
    /// Returns `None` until normalization succeeds or after it is invalidated.
    pub fn norm(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.normalization
                .as_ref()?
                .get_norm()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.normalization.as_ref()?.get_norm().cloned()
        }
    }

    /// Copy flattened normalized absorption on the current energy grid.
    /// Flattening removes the fitted post-edge trend; it is a display/analysis result,
    /// not the input used by AUTOBK. Returns `None` without valid normalization.
    pub fn flat(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.normalization
                .as_ref()?
                .get_flat()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.normalization.as_ref()?.get_flat().cloned()
        }
    }

    /// Copy the fitted pre-edge baseline on the current energy grid, in input mu units.
    /// Returns `None` without valid pre/post-edge normalization.
    pub fn pre_edge(&self) -> Option<DVector<f64>> {
        let normalization = self.normalization.as_ref()?;
        match normalization {
            normalization::NormalizationMethod::PrePostEdge(prepost) => {
                #[cfg(feature = "ndarray-compat")]
                {
                    prepost
                        .get_pre_edge()
                        .map(|x| DVector::from_vec(x.to_vec()))
                }
                #[cfg(not(feature = "ndarray-compat"))]
                {
                    prepost.get_pre_edge().cloned()
                }
            }
            _ => None,
        }
    }

    /// Copy the fitted post-edge normalization curve on the current energy grid,
    /// in input mu units. Returns `None` without valid pre/post-edge normalization.
    pub fn post_edge(&self) -> Option<DVector<f64>> {
        let normalization = self.normalization.as_ref()?;
        match normalization {
            normalization::NormalizationMethod::PrePostEdge(prepost) => {
                #[cfg(feature = "ndarray-compat")]
                {
                    prepost
                        .get_post_edge()
                        .map(|x| DVector::from_vec(x.to_vec()))
                }
                #[cfg(not(feature = "ndarray-compat"))]
                {
                    prepost.get_post_edge().cloned()
                }
            }
            _ => None,
        }
    }

    #[cfg(feature = "ndarray-compat")]
    /// Borrow the background wave-number grid in Å⁻¹ as an ndarray view.
    /// Available with `ndarray-compat`; returns `None` before a valid background result.
    pub fn k_view(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        self.background.as_ref()?.get_k_view()
    }

    #[cfg(feature = "ndarray-compat")]
    /// Borrow dimensionless, unweighted EXAFS as an ndarray view.
    /// Available with `ndarray-compat`; returns `None` before a valid background result.
    pub fn chi_view(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        self.background.as_ref()?.get_chi_view()
    }

    /// Borrow the currently stored forward-transform k-weight, if configured.
    /// Before `fft()` resolves its settings this can be an unresolved user value.
    pub fn kweight(&self) -> Option<&f64> {
        self.xftf.as_ref()?.get_kweight()
    }

    /// Calculate an owned `chi(k) * k^w` vector on the background grid.
    /// Here `w` is the currently stored forward-transform k-weight; units are Å⁻ʷ.
    /// This getter applies no window and does not resample onto the Larch FFT grid.
    /// Returns `None` if background results or the forward k-weight are unavailable.
    pub fn chi_kweighted(&self) -> Option<DVector<f64>> {
        let k = DVector::from_column_slice(self.k()?);
        let chi = DVector::from_column_slice(self.chi()?);
        let kweight = self.kweight()?;

        Some(chi.component_mul(&k.map(|x| x.powf(kweight.to_owned()))))
    }

    /// Borrow the complete stored one-sided complex real-FFT representation.
    /// Unlike the component getters, this includes frequencies beyond `rmax_out`.
    /// Returns `None` without a valid forward transform; no buffer is cloned.
    pub fn chir(&self) -> Option<&DynRealDft<f64>> {
        self.xftf.as_ref()?.get_chir()
    }

    /// Copy the forward-transform magnitude on [`Self::r`], limited by `rmax_out`.
    /// For dimensionless chi and k-weight w, units are Å⁻⁽ʷ⁺¹⁾. Returns `None` without
    /// a valid forward transform. These magnitudes are not normalized to a peak height.
    pub fn chir_mag(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftf
                .as_ref()?
                .get_chir_mag()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftf.as_ref()?.get_chir_mag().cloned()
        }
    }

    /// Wavenumbers for `kwin()`. In Larch mode these follow the resampled
    /// FFT grid, which can differ from the background's `k()` spacing.
    pub fn kwin_k(&self) -> Option<DVector<f64>> {
        let ft = self.xftf.as_ref()?;
        let len = ft.get_kwin()?.len();
        if ft.grid == super::xrayfft::FFTGrid::Larch {
            let step = *ft.get_kstep()?;
            Some(DVector::from_iterator(
                len,
                (0..len).map(|i| i as f64 * step),
            ))
        } else {
            Some(DVector::from_column_slice(self.k()?))
        }
    }

    /// Copy the dimensionless forward-transform window. Pair it with [`Self::kwin_k`],
    /// which can differ from the background grid in Larch mode. Returns `None` without
    /// a valid forward transform.
    pub fn kwin(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftf
                .as_ref()?
                .get_kwin()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftf.as_ref()?.get_kwin().cloned()
        }
    }

    /// Copy the real Fourier component on [`Self::r`], limited by `rmax_out`.
    /// Units are Å⁻⁽ʷ⁺¹⁾ for dimensionless chi and k-weight w; the forward exponent
    /// is negative. Returns `None` without a valid forward transform.
    pub fn chir_real(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftf
                .as_ref()?
                .get_chir_real()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftf.as_ref()?.get_chir_real()
        }
    }

    /// Copy the imaginary Fourier component on [`Self::r`], limited by `rmax_out`.
    /// Units are Å⁻⁽ʷ⁺¹⁾ for dimensionless chi and k-weight w; the forward exponent
    /// is negative. Returns `None` without a valid forward transform.
    pub fn chir_imag(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftf
                .as_ref()?
                .get_chir_imag()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftf.as_ref()?.get_chir_imag()
        }
    }

    /// Copy the reported Fourier distance grid in Å, limited by `rmax_out`.
    /// Scattering phase shifts are not corrected: a peak position is not directly a
    /// bond length. Returns `None` without a valid forward transform.
    pub fn r(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftf
                .as_ref()?
                .get_r()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftf.as_ref()?.get_r().cloned()
        }
    }

    /// Copy the inverse-transform wave-number grid in Å⁻¹, limited by `qmax_out`.
    /// Returns `None` without a valid inverse transform.
    pub fn q(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftr
                .as_ref()?
                .get_q()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftr.as_ref()?.get_q().cloned()
        }
    }

    /// Copy the real inverse-transform signal on [`Self::q`].
    /// Forward k-weighting and windowing remain in the signal; with inverse r-weight
    /// zero its units are Å⁻ʷ for forward weight w. Returns `None` without a valid
    /// inverse transform. This does not reconstruct removed background or lost data.
    pub fn chiq(&self) -> Option<DVector<f64>> {
        #[cfg(feature = "ndarray-compat")]
        {
            self.xftr
                .as_ref()?
                .get_chiq()
                .map(|x| DVector::from_vec(x.to_vec()))
        }
        #[cfg(not(feature = "ndarray-compat"))]
        {
            self.xftr.as_ref()?.get_chiq()
        }
    }
}

// Simple unit tests for this file.

#[cfg(test)]
pub mod tests {

    use super::*;
    use crate::xafs::io;
    use crate::xafs::tests::PARAM_LOADTXT;
    use crate::xafs::tests::TEST_TOL;
    use crate::xafs::tests::TEST_TOL_LESS_ACC;
    use crate::xafs::tests::TOP_DIR;
    use data_reader::reader::{load_txt_f64, Delimiter, ReaderParams};

    use approx::assert_abs_diff_eq;

    #[test]
    fn test_xafs_group_name_from_string() {
        let mut xafs_group = XASSpectrum::new();
        xafs_group.set_name("test".to_string());
        assert_eq!(xafs_group.name, Some("test".to_string()));
    }

    #[test]
    fn test_xafs_group_name_from_str() {
        let mut xafs_group = XASSpectrum::new();
        xafs_group.set_name("test");
        assert_eq!(xafs_group.name, Some("test".to_string()));

        let name = String::from("test");

        let mut xafs_group = XASSpectrum::new();
        xafs_group.set_name(name.clone());
        assert_eq!(xafs_group.name, Some("test".to_string()));

        println!("name: {}", name);
    }

    #[test]
    fn test_xafs_group_spectrum_from_vec() {
        let energy: Vec<f64> = vec![1.0, 2.0, 3.0];
        let mu: Vec<f64> = vec![4.0, 5.0, 6.0];
        let mut xafs_group = XASSpectrum::new();
        xafs_group.set_spectrum(energy, mu);
        assert_eq!(
            xafs_group.raw_energy,
            Some(DVector::from_vec(vec![1.0, 2.0, 3.0]))
        );
        assert_eq!(
            xafs_group.raw_mu,
            Some(DVector::from_vec(vec![4.0, 5.0, 6.0]))
        );
    }

    #[test]
    #[cfg(feature = "ndarray-compat")]
    fn test_xafs_group_normalization() {
        let test_file = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut xafs_group = io::load_spectrum_QAS_trans(&test_file).unwrap();

        let _ = xafs_group.normalize();

        let reference_path =
            String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS_pre_post_edge_expected.dat";
        let reference = load_txt_f64(&reference_path, &PARAM_LOADTXT).unwrap();

        let expected_norm = reference.get_col(4);

        let normalization = xafs_group.normalization.unwrap();
        let norm = normalization.get_norm().unwrap();
        norm.iter()
            .zip(expected_norm.iter())
            .for_each(|(x, y)| assert_abs_diff_eq!(x, y, epsilon = TEST_TOL_LESS_ACC));
    }

    #[test]
    #[cfg(not(feature = "ndarray-compat"))]
    fn test_xafs_group_normalization_nalgebra_smoke() {
        let test_file = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut xafs_group = io::load_spectrum_QAS_trans(&test_file).unwrap();

        xafs_group.normalize().unwrap();
        let norm = xafs_group
            .normalization
            .as_ref()
            .and_then(|method| method.get_norm())
            .unwrap();
        assert_eq!(norm.len(), xafs_group.energy.as_ref().unwrap().len());
        assert!(norm.iter().all(|value| value.is_finite()));
    }

    #[test]
    fn test_find_e0_rejects_non_monotonic_energy() {
        let mut spectrum = XASSpectrum::new();
        spectrum.energy = Some(DVector::from_vec(vec![1.0, 3.0, 2.0]));
        spectrum.mu = Some(DVector::from_vec(vec![1.0, 2.0, 3.0]));

        let err = spectrum.find_e0().unwrap_err();
        assert!(matches!(
            err,
            XAFSError::Data(DataError::NonMonotonicEnergy { .. })
        ));
    }

    #[test]
    fn test_normalize_rejects_non_finite_input() {
        let mut spectrum = XASSpectrum::new();
        spectrum.energy = Some(DVector::from_vec(vec![1.0, 2.0, 3.0]));
        spectrum.mu = Some(DVector::from_vec(vec![1.0, f64::NAN, 3.0]));

        let err = spectrum.normalize().unwrap_err();
        assert!(matches!(
            err,
            XAFSError::Data(DataError::NonFiniteValues { .. })
        ));
    }

    #[test]
    fn test_calc_background_rejects_length_mismatch() {
        let mut spectrum = XASSpectrum::new();
        spectrum.energy = Some(DVector::from_vec(vec![1.0, 2.0, 3.0]));
        spectrum.mu = Some(DVector::from_vec(vec![1.0, 2.0]));

        let err = spectrum.calc_background().unwrap_err();
        assert!(matches!(
            err,
            XAFSError::Data(DataError::LengthMismatch { .. })
        ));
    }

    #[test]
    fn test_interpolate_spectrum_updates_energy_and_mu() {
        let mut spectrum = XASSpectrum::new();
        spectrum.set_spectrum(vec![0.0, 1.0, 2.0, 3.0], vec![0.0, 2.0, 4.0, 6.0]);

        spectrum.interpolate_spectrum(vec![0.5, 1.5, 2.5]).unwrap();

        assert_eq!(
            spectrum.energy.as_ref().unwrap(),
            &DVector::from_vec(vec![0.5, 1.5, 2.5])
        );

        let mu = spectrum.mu.as_ref().unwrap();
        assert_abs_diff_eq!(mu[0], 1.0, epsilon = TEST_TOL);
        assert_abs_diff_eq!(mu[1], 3.0, epsilon = TEST_TOL);
        assert_abs_diff_eq!(mu[2], 5.0, epsilon = TEST_TOL);
    }

    #[test]
    fn test_interpolate_spectrum_missing_raw_mu_keeps_existing_mu() {
        let mut spectrum = XASSpectrum::new();
        spectrum.raw_energy = Some(DVector::from_vec(vec![0.0, 1.0]));
        spectrum.raw_mu = None;
        spectrum.mu = Some(DVector::from_vec(vec![42.0]));

        let err = spectrum.interpolate_spectrum(vec![0.25, 0.75]).unwrap_err();
        assert!(matches!(
            err,
            XAFSError::Data(DataError::MissingData { ref field }) if field == "raw_mu"
        ));

        assert_eq!(
            spectrum.energy.as_ref().unwrap(),
            &DVector::from_vec(vec![0.25, 0.75])
        );
        assert_eq!(
            spectrum.mu.as_ref().unwrap(),
            &DVector::from_vec(vec![42.0])
        );
    }

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

        let k_slice = spectrum.k().unwrap();
        let chi_slice = spectrum.chi().unwrap();
        let k_view = spectrum.k_view().unwrap();
        let chi_view = spectrum.chi_view().unwrap();

        assert_eq!(k_slice.len(), k_view.len());
        assert_eq!(chi_slice.len(), chi_view.len());

        for (slice, view) in k_slice.iter().zip(k_view.iter()) {
            assert_abs_diff_eq!(slice, view, epsilon = TEST_TOL);
        }
        for (slice, view) in chi_slice.iter().zip(chi_view.iter()) {
            assert_abs_diff_eq!(slice, view, epsilon = TEST_TOL_LESS_ACC);
        }

        spectrum.fft()?;
        assert!(spectrum.chir_mag().is_some());
        Ok(())
    }
}