rexafs 0.2.5

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

// Standard library dependencies
use std::cmp;
use std::error::Error;
// External dependencies
use nalgebra::DVector;
use ndarray::{Array, Array1, ArrayBase, ArrayView1, Axis, Ix1, OwnedRepr, Slice};
use serde::{Deserialize, Serialize};

// load dependencies
use super::bessel_i0;
use super::io;

// Load local traits
use super::mathutils::MathUtils;

// Load local functions
use crate::xafs::mathutils::index_of;

// Constants
/// Conventional small energy scale in eV used by energy-step utilities.
/// The default backend returns this fallback when an averaging slice is empty;
/// it is a numerical default, not a measured energy resolution.
pub const TINY_ENERGY: f64 = 0.005;

/// Physical constants used in rexafs
///
/// # Example
/// ```
/// use rexafs::xafs::xafsutils::constants;
///
/// assert_eq!(constants::h, 6.62607015e-34);
/// ```
#[path = "constants.rs"]
pub mod constants;

/// Convert excess photon energy E−E0 (eV) and photoelectron wave number k (Å⁻¹).
///
/// The nonrelativistic relation is E−E0 = hbar²*k²/(2*m_e), with unit conversion
/// encoded by `constants::KTOE` ≈ 3.809982110968585 eV Ų. Subtract E0 before
/// calling `etok`; these helpers do not know the absorption-edge energy.
/// Vector implementations allocate results and leave borrowed inputs unchanged.
pub trait XAFSUtils {
    /// Return sqrt((E−E0)*ETOK) in Å⁻¹. Scalar, Vec and DVector values below
    /// zero are clamped to zero; NaNs propagate. The legacy ndarray-array
    /// implementation instead takes the raw square root and returns NaN for
    /// negative excess energy. AUTOBK uses a separate signed interpolation grid.
    fn etok(&self) -> Self;
    /// Return k²*KTOE in eV; E0 is not added. Negative k is squared, so its
    /// sign is lost. No finiteness validation is performed.
    fn ktoe(&self) -> Self;
}

impl XAFSUtils for f64 {
    fn etok(&self) -> Self {
        if *self < 0.0 {
            return 0.0;
        }

        (self * constants::ETOK).sqrt()
    }

    fn ktoe(&self) -> Self {
        self.powi(2) * constants::KTOE
    }
}

impl XAFSUtils for Vec<f64> {
    fn etok(&self) -> Self {
        self.iter().map(|x| x.etok()).collect()
    }

    fn ktoe(&self) -> Self {
        self.iter().map(|x| x.ktoe()).collect()
    }
}

impl XAFSUtils for nalgebra::DVector<f64> {
    fn etok(&self) -> Self {
        self.map(|x| {
            if x < 0.0 {
                0.0
            } else {
                (x * constants::ETOK).sqrt()
            }
        })
    }

    fn ktoe(&self) -> Self {
        self.map(|x| x.powi(2) * constants::KTOE)
    }
}

impl XAFSUtils for ArrayBase<OwnedRepr<f64>, Ix1> {
    fn etok(&self) -> Self {
        self.mapv(|x| (x * constants::ETOK).sqrt())
    }

    fn ktoe(&self) -> Self {
        self.mapv(|x| x.powi(2) * constants::KTOE)
    }
}

/// Line-profile family used for smoothing by convolution.
/// Widths use the same units as the supplied x axis; these kernels are not the
/// Fourier-window families in `FTWindow`. See
/// [SciPy's Voigt-profile definition](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.voigt_profile.html).
#[derive(Debug, Clone, Copy, Default)]
pub enum ConvolveForm {
    /// Lorentzian kernel; sigma is its half width at half maximum.
    #[default]
    Lorentzian,
    /// Gaussian kernel; sigma is its standard deviation.
    Gaussian,
    /// Gaussian–Lorentzian convolution; sigma is the Gaussian standard deviation
    /// and gamma is the Lorentzian half width at half maximum.
    Voigt,
}
/// Smooth y(x) by interpolation, reflected extension and direct convolution.
///
/// Returns a newly allocated array on the original x coordinates, in the same
/// y units. Widths sigma (default 1.0) and gamma (default sigma) use x units.
/// Lorentzian uses sigma as half width at half maximum and ignores gamma;
/// Gaussian uses sigma as standard deviation; Voigt uses both widths.
/// Increasing the width suppresses narrow features as well as noise.
///
/// The requested xstep defaults to the smallest successive x difference.
/// npad defaults to 5 and extends the interpolation domain by that many requested
/// steps on each side. The code reflects the interpolated signal, constructs a
/// sampled kernel, normalizes its sum, and evaluates a direct valid convolution;
/// it does not call an FFT. See the convolution definition in
/// [NumPy's reference](https://numpy.org/doc/stable/reference/generated/numpy.convolve.html).
/// The interpolation array is capped at 50 times the input length. If that cap
/// is reached, its actual spacing differs from xstep while kernel widths still
/// use the requested step; avoid excessively fine xstep values.
///
/// Use matching finite arrays on a strictly increasing x axis and positive finite
/// widths and spacing. Interpolation or invalid convolution dimensions can return
/// an error. This low-level helper does not validate every nonfinite value.
/// The ndarray compatibility helper takes ownership through `Into<Array1>`.
/// It does not clamp npad or validate short/mismatched arrays before indexing;
/// short arrays, invalid ranges or xstep below 1e-12 can panic.
pub fn smooth<T: Into<Array1<f64>>>(
    x: T,
    y: T,
    sigma: Option<f64>,
    gamma: Option<f64>,
    xstep: Option<f64>,
    npad: Option<i32>,
    conv_form: ConvolveForm,
) -> Result<Array1<f64>, Box<dyn Error>> {
    const TINY: f64 = 1e-12;

    let x: Array1<f64> = x.into();
    let y: Array1<f64> = y.into();
    let npad = npad.unwrap_or(5);

    let x_diff = x.diff();
    let xstep = xstep.unwrap_or(x_diff.min());

    if xstep < TINY {
        todo!("Cannot smooth data: must be strictly increasing. Impliment error handling");
    }

    let sigma = sigma.unwrap_or(1.0);
    let gamma = gamma.unwrap_or(sigma);

    let xmin = xstep * ((x.min() - npad as f64 * xstep) / xstep).floor();
    let xmax = xstep * ((x.max() + npad as f64 * xstep) / xstep).floor();
    let npts1 = 1 + ((xmax - xmin + xstep * 0.1) / xstep).abs() as i32;
    let npts = npts1.min(50 * x.len() as i32);

    let x0: Array1<f64> = Array1::linspace(xmin, xmax, npts as usize);
    let y0: Array1<f64> = if let (Some(x_slice), Some(y_slice)) = (x.as_slice(), y.as_slice()) {
        x0.interpolate(x_slice, y_slice)?
    } else {
        x0.interpolate(&x.to_vec(), &y.to_vec())?
    };

    let sigma = sigma / xstep;
    let gamma = gamma / xstep;

    let wx: Array1<f64> = Array1::range(0.0, 2.0 * npts as f64, 1.0);
    let win: Array1<f64> = match conv_form {
        ConvolveForm::Gaussian => wx.gaussian(npts as f64, sigma),
        ConvolveForm::Voigt => wx.voigt(npts as f64, sigma, gamma),
        ConvolveForm::Lorentzian => wx.lorentzian(npts as f64, sigma),
    };

    let y1 = ndarray::concatenate(
        ndarray::Axis(0),
        &[
            y0.slice_axis(Axis(0), Slice::from(0..npts).step_by(-1)),
            y0.view(),
            y0.slice_axis(Axis(0), Slice::from((-npts as i32)..-1).step_by(-1)),
        ],
    )?;

    let y2 = convolve_valid(&y1, &(&win / win.sum()))?;

    let y2 = if y2.len() > x0.len() {
        let nex = ((y2.len() - x0.len()) / 2) as usize;
        let y2 = y2.slice_axis(Axis(0), Slice::from(nex..(nex + x0.len())).step_by(1));
        y2
    } else {
        y2.view()
    };

    Ok(
        if let (Some(x0_slice), Some(y2_slice)) = (x0.as_slice(), y2.as_slice()) {
            x.interpolate(x0_slice, y2_slice)?
        } else {
            x.interpolate(&x0.to_vec(), &y2.to_vec())?
        },
    )
}

fn convolve_valid(
    signal: &Array1<f64>,
    kernel: &Array1<f64>,
) -> Result<Array1<f64>, Box<dyn Error>> {
    if signal.is_empty() || kernel.is_empty() || signal.len() < kernel.len() {
        return Err(std::io::Error::new(
            std::io::ErrorKind::InvalidInput,
            "invalid convolution input lengths",
        )
        .into());
    }

    let out_len = signal.len() - kernel.len() + 1;
    let mut out = Array1::zeros(out_len);

    for i in 0..out_len {
        let mut acc = 0.0;
        for j in 0..kernel.len() {
            // Match convolution semantics used by fft-based implementation.
            acc += signal[i + j] * kernel[kernel.len() - 1 - j];
        }
        out[i] = acc;
    }

    Ok(out)
}

/// Nudge nearly repeated coordinates without deleting any samples.
///
/// Takes ownership through `Into<Array1>` and returns an array with the same length. For a difference d
/// between the current and preceding non-NaN original values, if abs(d) < tiny,
/// the added offset accumulates by max(tiny, frac * abs(d)). Defaults are
/// tiny=1e-7 in coordinate units and dimensionless frac=1e-6. The increment does
/// not use the next sample. This is a numerical coordinate adjustment, not an
/// averaging or merging operation on paired absorption values.
///
/// `sort=Some(true)` sorts a copy first; None/false retain the input order.
/// Use finite values when sorting: the floating-point comparison unwraps and
/// panics on NaN. Without sorting, NaNs are retained and skipped when tracking
/// the previous value. This helper neither guarantees strict monotonicity after
/// nudging nor changes an associated mu array.
#[cfg(feature = "ndarray-compat")]
pub fn remove_dups_array1<T: Into<ArrayBase<OwnedRepr<f64>, Ix1>>>(
    arr: T,
    tiny: Option<f64>,
    frac: Option<f64>,
    sort: Option<bool>,
) -> ArrayBase<OwnedRepr<f64>, Ix1> {
    let mut arr = arr.into();
    let tiny = tiny.unwrap_or(1e-7);
    let frac = frac.unwrap_or(1e-6);

    if arr.len() < 2 {
        return arr;
    }

    if let Some(true) = sort {
        let mut arr_sort = arr.to_vec();
        arr_sort.sort_by(|a, b| a.partial_cmp(b).unwrap());
        arr = Array1::from_vec(arr_sort);
    }

    let mut previous_value = f64::NAN;
    let mut previous_add = 0.0;

    let mut add = Array1::zeros(arr.len());

    for i in 1..arr.len() {
        if !arr[i - 1].is_nan() {
            previous_value = arr[i - 1];
            previous_add = add[i - 1];
        }
        let value = arr[i];
        if value.is_nan() || previous_value.is_nan() {
            continue;
        }
        let diff = (value - previous_value).abs();
        if diff < tiny {
            add[i] = previous_add + f64::max(tiny, frac * diff);
        }
    }

    arr = arr + add;

    arr
}

/// Nudge nearly repeated coordinates without deleting any samples.
///
/// Returns a newly allocated vector with the same length. For a difference d
/// between the current and preceding non-NaN original values, if abs(d) < tiny,
/// the added offset accumulates by max(tiny, frac * abs(d)). Defaults are
/// tiny=1e-7 in coordinate units and dimensionless frac=1e-6. The increment does
/// not use the next sample. This is a numerical coordinate adjustment, not an
/// averaging or merging operation on paired absorption values.
///
/// `sort=Some(true)` sorts a copy first; None/false retain the input order.
/// Use finite values when sorting: the floating-point comparison unwraps and
/// panics on NaN. Without sorting, NaNs are retained and skipped when tracking
/// the previous value. This helper neither guarantees strict monotonicity after
/// nudging nor changes an associated mu array.
pub fn remove_dups(
    arr: &DVector<f64>,
    tiny: Option<f64>,
    frac: Option<f64>,
    sort: Option<bool>,
) -> DVector<f64> {
    let tiny = tiny.unwrap_or(1e-7);
    let frac = frac.unwrap_or(1e-6);

    if arr.len() < 2 {
        return arr.clone();
    }

    let arr = if let Some(true) = sort {
        let mut arr_sort = arr.as_slice().to_vec();
        arr_sort.sort_by(|a, b| a.partial_cmp(b).unwrap());
        DVector::from_vec(arr_sort)
    } else {
        arr.clone()
    };

    let mut previous_value = f64::NAN;
    let mut previous_add = 0.0;

    let mut add = DVector::zeros(arr.len());

    for i in 1..arr.len() {
        if !arr[i - 1].is_nan() {
            previous_value = arr[i - 1];
            previous_add = add[i - 1];
        }
        let value = arr[i];
        if value.is_nan() || previous_value.is_nan() {
            continue;
        }
        let diff = (value - previous_value).abs();
        if diff < tiny {
            add[i] = previous_add + f64::max(tiny, frac * diff);
        }
    }

    arr + add
}

/// Allocate paired arrays retaining only entries where both values are finite.
///
/// Despite the historical name, infinities are removed as well as NaNs. The
/// inputs are borrowed views; iteration stops at the shorter input, so unequal
/// lengths silently discard an unmatched tail. Units and retained order are unchanged.
pub fn remove_nan2(
    arr1: ArrayView1<'_, f64>,
    arr2: ArrayView1<'_, f64>,
) -> (Array1<f64>, Array1<f64>) {
    let (arr1, arr2): (Vec<f64>, Vec<f64>) = arr1
        .iter()
        .zip(arr2.iter())
        .filter(|(e, m)| e.is_finite() && m.is_finite())
        .unzip();

    (arr1.into(), arr2.into())
}

/// Estimate a small representative energy interval in eV.
///
/// Defaults are frac_ignore=0.01, nave=10 and sort=false. Differences between
/// successive energy values are sorted. With n input energies, start is
/// floor(frac_ignore*n), and end is min(start+nave, n-2). The mean uses the
/// half-open difference slice [start,end), so the largest interval is excluded.
/// The ignored count is based on energy sample count, not difference count.
///
/// This heuristic supplies smoothing scales for edge finding; it is not an
/// energy calibration or an uncertainty estimate. Borrowed input is unchanged.
/// Use finite energies and finite nonnegative fraction settings: comparison of
/// NaNs can panic. Sorting, when requested, operates on a copy.
/// The input is consumed through `Into<Array1>`. This legacy helper does not use
/// the default backend's 0.005 eV fallback. A short array or unusable averaging
/// slice can panic or return NaN.
#[cfg(feature = "ndarray-compat")]
pub fn find_energy_step_array1<T: Into<ArrayBase<OwnedRepr<f64>, Ix1>>>(
    energy: T,
    frac_ignore: Option<f64>,
    nave: Option<usize>,
    sort: Option<bool>,
) -> f64 {
    let mut energy = energy.into();

    if let Some(true) = sort {
        let mut energy_sort = energy.to_vec();
        energy_sort.sort_by(|a, b| a.partial_cmp(b).unwrap());
        energy = Array1::from_vec(energy_sort);
    }

    let frac_ignore = frac_ignore.unwrap_or(0.01);
    let nave = nave.unwrap_or(10);
    let mut ediff = (&energy.slice(ndarray::s![1..]) - &energy.slice(ndarray::s![..-1]))
        .to_owned()
        .to_vec();

    let nskip = (frac_ignore * energy.len() as f64) as usize;

    ediff.sort_by(|a, b| a.partial_cmp(b).unwrap());

    let ediff_end = cmp::min(nskip + nave, ediff.len() - 1);

    return ediff[nskip..ediff_end].iter().sum::<f64>() / (ediff_end - nskip) as f64;
}
/// Estimate an absorption edge in eV from the derivative of mu(E).
///
/// Runs `_find_e0` on the full spectrum, then refines a neighborhood extending
/// up to 75 samples on each side of the candidate. A peak requires adjacent
/// high-derivative samples to reduce isolated-glitch sensitivity. This is a
/// numerical edge estimate, not independent energy calibration or evidence that
/// an edge is physically unique. Inspect noisy, multi-edge or narrow scans.
/// See [Larch's edge-finding reference](https://xraypy.github.io/xraylarch/xafs_preedge.html#the-find-e0-function)
/// for the method's purpose; detailed thresholds and fallback rules here are
/// implementation choices.
/// This legacy ndarray routine consumes its arrays and performs a smoothed
/// second pass. The numerical result can differ from the default backend,
/// including endpoint handling; no universal Larch equivalence is promised.
/// Unlike the default backend, invalid or short refinement ranges can panic.
#[cfg(feature = "ndarray-compat")]
pub fn find_e0_array1<T: Into<ArrayBase<OwnedRepr<f64>, Ix1>>>(
    energy: T,
    mu: T,
) -> Result<f64, Box<dyn Error>> {
    let energy: ArrayBase<OwnedRepr<f64>, Ix1> = energy.into();
    let mu: ArrayBase<OwnedRepr<f64>, Ix1> = mu.into();

    let (e1, ie0, estep) = _find_e0_array1(energy.clone(), mu.clone(), None, None)?;
    let istart = (ie0 as i32 - 75).max(2) as usize;
    let istop = (ie0 + 75).min(energy.len() - 2);

    let (mut e0, ix, ex) = _find_e0_array1(
        energy.slice(ndarray::s![istart..istop]).to_owned(),
        mu.slice(ndarray::s![istart..istop]).to_owned(),
        Some(estep),
        Some(true),
    )?;

    if ix < 1 {
        e0 = energy[istart + 2];
    }

    Ok(e0)
}

/// Find one derivative-based edge candidate and return (energy_eV, index, step_eV).
///
/// The estimate uses gradient(mu)/gradient(energy) after nudging duplicate
/// energy coordinates. A normalized derivative threshold starts at 0.60 for
/// more than 20 samples and 0.30 otherwise, and can be halved twice. A selected
/// peak must also have neighboring samples above the threshold; end regions
/// are excluded. These numerical thresholds are rexafs choices, not confidence
/// levels. If no candidate passes, the initial index zero can be returned.
///
/// `estep=None` uses half `find_energy_step`; `use_smooth=None` is false.
/// Smoothing, when used, is Lorentzian with width 3*estep and sample spacing
/// estep. The index refers to the supplied sample order; duplicate nudging may
/// shift the returned energy slightly. Finite increasing energy in eV and
/// matched absorption samples are the intended inputs.
/// This historical ndarray implementation consumes its arrays and lacks
/// comprehensive length/finiteness/range checks. Optional smoothing is performed,
/// but its error is unwrapped; malformed input can therefore panic.
/// The DVector helper under ndarray-compat does not perform that smoothing.
#[cfg(feature = "ndarray-compat")]
pub fn _find_e0_array1<T: Into<ArrayBase<OwnedRepr<f64>, Ix1>> + Clone>(
    energy: T,
    mu: T,
    estep: Option<f64>,
    use_smooth: Option<bool>,
) -> Result<(f64, usize, f64), Box<dyn Error>> {
    let en: ArrayBase<OwnedRepr<f64>, Ix1> =
        remove_dups_array1(energy.clone().into(), None, None, None);
    let mu: ArrayBase<OwnedRepr<f64>, Ix1> = mu.into();

    let estep =
        estep.unwrap_or(find_energy_step_array1(energy.clone(), None, None, Some(false)) / 2.0);

    let nmin = 2.max(en.len() / 100);

    let dmu: ArrayBase<OwnedRepr<f64>, Ix1> = if let Some(true) = use_smooth {
        // todo!("smooth not implemented yet");
        smooth(
            energy.into(),
            mu.gradient() / en.gradient(),
            Some(3.0 * estep),
            None,
            Some(estep),
            None,
            ConvolveForm::Lorentzian,
        )
        .unwrap()
    } else {
        mu.gradient() / en.gradient()
    };

    let dmin = *dmu
        .slice(ndarray::s![(nmin as i32)..(1 - nmin as i32)])
        .iter()
        .map(|a| if a.is_finite() { a } else { &-1.0 })
        .min_by(|a, b| a.partial_cmp(b).unwrap())
        .unwrap();

    let dm_ptp = dmu
        .slice(ndarray::s![(nmin as i32)..(1 - nmin as i32)])
        .to_vec()
        .ptp();

    let dmu = (dmu - dmin) / dm_ptp;

    let mut dhigh = if en.len() > 20 { 0.60 } else { 0.30 };

    let mut high_deriv_pts: Vec<usize> = dmu
        .indexed_iter()
        .filter(|(_, a)| a > &&dhigh)
        .map(|(i, _)| i)
        .collect();

    if high_deriv_pts.len() < 3 {
        for _ in 0..2 {
            if high_deriv_pts.len() > 3 {
                break;
            }

            dhigh *= 0.5;

            high_deriv_pts = dmu
                .indexed_iter()
                .filter(|(_, a)| a > &&dhigh)
                .map(|(i, _)| i)
                .collect();
        }
    }

    if high_deriv_pts.len() < 3 {
        high_deriv_pts = dmu
            .indexed_iter()
            .filter(|(_, a)| a.is_finite())
            .take(1)
            .map(|(i, _)| i)
            .collect();
    }

    let mut high_deriv_mask = vec![false; dmu.len()];
    for &idx in &high_deriv_pts {
        if idx < high_deriv_mask.len() {
            high_deriv_mask[idx] = true;
        }
    }

    let mut imax = 0;
    let mut dmax = 0.0;

    for i in &high_deriv_pts {
        if i < &nmin || i > &(dmu.len() - nmin) {
            continue;
        }

        let idx = *i;
        let has_prev = idx > 0 && high_deriv_mask[idx - 1];
        let has_next = idx + 1 < high_deriv_mask.len() && high_deriv_mask[idx + 1];
        if dmu[idx] > dmax && has_prev && has_next {
            dmax = dmu[idx];
            imax = *i;
        }
    }

    Ok((en[imax], imax, estep))
}

/// Estimate a small representative energy interval in eV.
///
/// Defaults are frac_ignore=0.01, nave=10 and sort=false. Differences between
/// successive energy values are sorted. With n input energies, start is
/// floor(frac_ignore*n), and end is min(start+nave, n-2). The mean uses the
/// half-open difference slice [start,end), so the largest interval is excluded.
/// The ignored count is based on energy sample count, not difference count.
///
/// This heuristic supplies smoothing scales for edge finding; it is not an
/// energy calibration or an uncertainty estimate. Borrowed input is unchanged.
/// Use finite energies and finite nonnegative fraction settings: comparison of
/// NaNs can panic. Sorting, when requested, operates on a copy.
/// This compatibility helper does not supply the default backend's empty-slice
/// fallback. Empty/short arrays or an invalid averaging slice can panic or
/// return NaN; provide a nonempty usable averaging slice.
pub fn find_energy_step(
    energy: &DVector<f64>,
    frac_ignore: Option<f64>,
    nave: Option<usize>,
    sort: Option<bool>,
) -> f64 {
    let energy = if let Some(true) = sort {
        let mut energy_sort = energy.as_slice().to_vec();
        energy_sort.sort_by(|a, b| a.partial_cmp(b).unwrap());
        DVector::from_vec(energy_sort)
    } else {
        energy.clone()
    };

    let frac_ignore = frac_ignore.unwrap_or(0.01);
    let nave = nave.unwrap_or(10);

    // Calculate differences: energy[1..] - energy[..-1]
    let mut ediff: Vec<f64> = Vec::with_capacity(energy.len() - 1);
    for i in 1..energy.len() {
        ediff.push(energy[i] - energy[i - 1]);
    }

    let nskip = (frac_ignore * energy.len() as f64) as usize;

    ediff.sort_by(|a, b| a.partial_cmp(b).unwrap());

    let ediff_end = cmp::min(nskip + nave, ediff.len() - 1);

    ediff[nskip..ediff_end].iter().sum::<f64>() / (ediff_end - nskip) as f64
}

/// Find one derivative-based edge candidate and return (energy_eV, index, step_eV).
///
/// The estimate uses gradient(mu)/gradient(energy) after nudging duplicate
/// energy coordinates. A normalized derivative threshold starts at 0.60 for
/// more than 20 samples and 0.30 otherwise, and can be halved twice. A selected
/// peak must also have neighboring samples above the threshold; end regions
/// are excluded. These numerical thresholds are rexafs choices, not confidence
/// levels. If no candidate passes, the initial index zero can be returned.
///
/// `estep=None` uses half `find_energy_step`; `use_smooth=None` is false.
/// Smoothing, when used, is Lorentzian with width 3*estep and sample spacing
/// estep. The index refers to the supplied sample order; duplicate nudging may
/// shift the returned energy slightly. Finite increasing energy in eV and
/// matched absorption samples are the intended inputs.
/// The ndarray-feature DVector variant currently ignores use_smooth and always
/// uses the unsmoothed derivative. Its short-array/range checks are historical
/// and can panic; it is not the same routine as `_find_e0_array1`.
pub fn _find_e0(
    energy: &DVector<f64>,
    mu: &DVector<f64>,
    estep: Option<f64>,
    use_smooth: Option<bool>,
) -> Result<(f64, usize, f64), Box<dyn Error>> {
    let en = remove_dups(energy, None, None, None);

    let estep = estep.unwrap_or(find_energy_step(energy, None, None, Some(false)) / 2.0);

    let nmin = 2.max(en.len() / 100);

    // Calculate gradient: mu.gradient() / en.gradient()
    let mu_grad = mu.gradient();
    let en_grad = en.gradient();

    let dmu: DVector<f64> = if let Some(true) = use_smooth {
        // For now, skip smooth implementation - would need DVector version
        // Use simple gradient instead
        DVector::from_fn(mu_grad.len(), |i, _| {
            if en_grad[i].abs() > 1e-12 {
                mu_grad[i] / en_grad[i]
            } else {
                0.0
            }
        })
    } else {
        DVector::from_fn(mu_grad.len(), |i, _| {
            if en_grad[i].abs() > 1e-12 {
                mu_grad[i] / en_grad[i]
            } else {
                0.0
            }
        })
    };

    // Calculate dmin from the middle section
    let dmin = *dmu
        .as_slice()
        .iter()
        .skip(nmin)
        .take(dmu.len() - 2 * nmin)
        .filter(|a| a.is_finite())
        .min_by(|a, b| a.partial_cmp(b).unwrap())
        .unwrap_or(&-1.0);

    // Calculate peak-to-peak for normalization
    let middle_slice: Vec<f64> = dmu
        .as_slice()
        .iter()
        .skip(nmin)
        .take(dmu.len() - 2 * nmin)
        .cloned()
        .collect();
    let dm_ptp = middle_slice.ptp();

    // Normalize dmu
    let dmu = DVector::from_fn(dmu.len(), |i, _| (dmu[i] - dmin) / dm_ptp);

    let mut dhigh = if en.len() > 20 { 0.60 } else { 0.30 };

    let mut high_deriv_pts: Vec<usize> = dmu
        .iter()
        .enumerate()
        .filter(|(_, a)| a > &&dhigh)
        .map(|(i, _)| i)
        .collect();

    if high_deriv_pts.len() < 3 {
        for _ in 0..2 {
            if high_deriv_pts.len() > 3 {
                break;
            }

            dhigh *= 0.5;

            high_deriv_pts = dmu
                .iter()
                .enumerate()
                .filter(|(_, a)| a > &&dhigh)
                .map(|(i, _)| i)
                .collect();
        }
    }

    if high_deriv_pts.len() < 3 {
        high_deriv_pts = dmu
            .iter()
            .enumerate()
            .filter(|(_, a)| a.is_finite())
            .take(1)
            .map(|(i, _)| i)
            .collect();
    }

    let mut high_deriv_mask = vec![false; dmu.len()];
    for &idx in &high_deriv_pts {
        if idx < high_deriv_mask.len() {
            high_deriv_mask[idx] = true;
        }
    }

    let mut imax = 0;
    let mut dmax = 0.0;

    for i in &high_deriv_pts {
        if i < &nmin || i > &(dmu.len() - nmin) {
            continue;
        }

        let idx = *i;
        let has_prev = idx > 0 && high_deriv_mask[idx - 1];
        let has_next = idx + 1 < high_deriv_mask.len() && high_deriv_mask[idx + 1];
        if dmu[idx] > dmax && has_prev && has_next {
            dmax = dmu[idx];
            imax = *i;
        }
    }

    Ok((en[imax], imax, estep))
}

/// Estimate an absorption edge in eV from the derivative of mu(E).
///
/// Runs `_find_e0` on the full spectrum, then refines a neighborhood extending
/// up to 75 samples on each side of the candidate. A peak requires adjacent
/// high-derivative samples to reduce isolated-glitch sensitivity. This is a
/// numerical edge estimate, not independent energy calibration or evidence that
/// an edge is physically unique. Inspect noisy, multi-edge or narrow scans.
/// See [Larch's edge-finding reference](https://xraypy.github.io/xraylarch/xafs_preedge.html#the-find-e0-function)
/// for the method's purpose; detailed thresholds and fallback rules here are
/// implementation choices.
/// This ndarray-feature DVector entry point uses its unsmoothed legacy helper
/// in both passes. It lacks the default backend's short-window fallback and
/// can panic on short or malformed arrays. Inputs are borrowed unchanged.
pub fn find_e0(energy: &DVector<f64>, mu: &DVector<f64>) -> Result<f64, Box<dyn Error>> {
    let (e1, ie0, estep) = _find_e0(energy, mu, None, None)?;
    let istart = (ie0 as i32 - 75).max(2) as usize;
    let istop = (ie0 + 75).min(energy.len() - 2);

    // Extract slice as new DVector
    let energy_slice = DVector::from_iterator(
        istop - istart,
        energy.as_slice()[istart..istop].iter().cloned(),
    );
    let mu_slice =
        DVector::from_iterator(istop - istart, mu.as_slice()[istart..istop].iter().cloned());

    let (mut e0, ix, ex) = _find_e0(&energy_slice, &mu_slice, Some(estep), Some(true))?;

    if ix < 1 {
        e0 = energy[istart + 2];
    }

    Ok(e0)
}

/// Dimensionless Fourier-window families.
///
/// Windows reduce ringing at a truncated data range but also change amplitudes
/// and effective resolution. Their parameters are not interchangeable; see
/// [`ftwindow`] and the
/// [Larch window reference](https://xraypy.github.io/xraylarch/xafs_fourier.html#ftwindow-generating-fourier-transform-windows).
#[derive(Debug, Clone, Copy, Default, PartialEq, Serialize, Deserialize)]
pub enum FTWindow {
    #[default]
    /// Cosine-squared low/high tapers with an interior plateau; the helper default.
    Hanning,
    /// Linear low/high tapers with an interior plateau.
    Parzen,
    /// Quadratic low/high tapers with an interior plateau.
    Welch,
    /// Gaussian centered on the selected domain; `dx` is its standard deviation
    /// in axis units. Tails are evaluated over the complete supplied grid.
    Gaussian,
    /// A single sine arch across the selected domain.
    Sine,
    /// Modified Bessel window. `dx` controls both the domain and the numerical
    /// shape parameter; this is Larch's modified Kaiser convention, not unit-area scaling.
    KaiserBessel,
    /// Cosine-squared window with fractional taper geometry. Use dimensionless
    /// fractions for `dx`/`dx2`; equal numeric widths differ from ordinary Hanning.
    FHanning,
}

impl FTWindow {
    /// Construct this window on x; see `ftwindow` for units, defaults and grid assumptions.
    pub fn window(
        &self,
        x: &ArrayBase<OwnedRepr<f64>, Ix1>,
        xmin: Option<f64>,
        xmax: Option<f64>,
        dx: Option<f64>,
        dx2: Option<f64>,
    ) -> Result<Array1<f64>, Box<dyn Error>> {
        ftwindow(x, xmin, xmax, dx, dx2, Some(*self))
    }
}

/// Construct a dimensionless Fourier window on the supplied axis.
///
/// For k-space, x is in Å⁻¹; for R-space, x is in Å. Missing boundaries use
/// the axis extrema, `dx` defaults to 1, `dx2` follows `dx`, and a missing family
/// uses Hanning. Ordinary taper parameters have axis units. Gaussian uses `dx`
/// as its standard deviation, Kaiser–Bessel also uses its numeric value as a
/// shape parameter, and fractional Hanning uses fractions for taper geometry.
/// No division by area or by the sum of the samples is performed.
///
/// This low-level routine follows Larch's index geometry, intended for a
/// zero-origin uniform increasing axis and valid boundaries/tapers inside its
/// indexable range. It does not validate arbitrary grids or boundary settings;
/// invalid indices can panic. Prefer checked `XrayFFTF::xftf`/`XrayFFTR::xftr`
/// settings for the spectrum pipeline. This legacy helper requires at least
/// two samples and does not handle empty/single-sample grids. Input is unchanged.
///
/// The formulas and endpoint rounding follow the
/// [pinned Larch implementation](https://github.com/xraypy/xraylarch/blob/860d8a690c81eefb0e61dee4ca3703ef4b67e93d/larch/xafs/xafsft.py#L41).
pub fn ftwindow(
    x: &ArrayBase<OwnedRepr<f64>, Ix1>,
    xmin: Option<f64>,
    xmax: Option<f64>,
    dx: Option<f64>,
    dx2: Option<f64>,
    window: Option<FTWindow>,
) -> Result<Array1<f64>, Box<dyn Error>> {
    let window = window.unwrap_or_default();

    let mut dx1 = dx.unwrap_or(1.0);
    let mut dx2 = dx2.unwrap_or(dx1);

    let xmin = xmin.unwrap_or(x.min());
    let xmax = xmax.unwrap_or(x.max());

    let xstep = (x[x.len() - 1] - x[0]) / (x.len() as f64 - 1.0);
    let xeps = &xstep * 1e-4;

    let mut x1 = x.min().max(xmin - dx1 / 2.0);
    let mut x2 = xmin + dx1 / 2.0 + xeps;
    let mut x3 = xmax - dx2 / 2.0 - xeps;
    let mut x4 = x.max().min(xmax + dx2 / 2.0);

    let asint = |val: &f64| ((val + xeps) / xstep) as i32;

    match window {
        FTWindow::Gaussian => {
            dx1 = dx1.max(xeps);
        }

        FTWindow::FHanning => {
            if dx1 < 0.0 {
                dx1 = 0.0;
            }
            if dx2 > 1.0 {
                dx2 = 1.0;
            }
            x2 = x1 + xeps + dx1 * (xmax - xmin) / 2.0;
            x3 = x4 - xeps - dx2 * (xmax - xmin) / 2.0;
        }
        _ => {}
    }

    let (mut i1, mut i2, mut i3, mut i4) = (asint(&x1), asint(&x2), asint(&x3), asint(&x4));
    i1 = i1.max(0);
    i2 = i2.max(0);
    i3 = i3.min((x.len() - 1) as i32);
    i4 = i4.min((x.len() - 1) as i32);

    if i1 == i2 {
        i1 = (i2 - 1).max(0);
    }

    if i3 == i4 {
        i3 = (i4 - 1).max(i2);
    }

    (x1, x2, x3, x4) = (
        x[i1 as usize],
        x[i2 as usize],
        x[i3 as usize],
        x[i4 as usize],
    );
    if x1 == x2 {
        x2 += xeps;
    }

    if x3 == x4 {
        x4 += xeps;
    }

    let mut fwin = Array1::zeros(x.len());

    if i3 > i2 {
        fwin.slice_mut(ndarray::s![i2..i3]).fill(1.0);
    }

    match window {
        FTWindow::Hanning | FTWindow::FHanning => {
            fwin.slice_mut(ndarray::s![i1..=i2])
                .assign(&x.slice(ndarray::s![i1..=i2]).mapv(|x| {
                    (std::f64::consts::PI / 2.0 * (x - x1) / (x2 - x1))
                        .sin()
                        .powi(2)
                }));
            fwin.slice_mut(ndarray::s![i3..=i4])
                .assign(&x.slice(ndarray::s![i3..=i4]).mapv(|x| {
                    (std::f64::consts::PI / 2.0 * (x - x3) / (x4 - x3))
                        .cos()
                        .powi(2)
                }));
        }
        FTWindow::Parzen => {
            fwin.slice_mut(ndarray::s![i1..=i2])
                .assign(&x.slice(ndarray::s![i1..=i2]).mapv(|x| (x - x1) / (x2 - x1)));
            fwin.slice_mut(ndarray::s![i3..=i4]).assign(
                &x.slice(ndarray::s![i3..=i4])
                    .mapv(|x| 1.0 - (x - x3) / (x4 - x3)),
            );
        }
        FTWindow::Welch => {
            fwin.slice_mut(ndarray::s![i1..=i2]).assign(
                &x.slice(ndarray::s![i1..=i2])
                    .mapv(|x| 1.0 - ((x - x2) / (x2 - x1)).powi(2)),
            );
            fwin.slice_mut(ndarray::s![i3..=i4]).assign(
                &x.slice(ndarray::s![i3..=i4])
                    .mapv(|x| 1.0 - ((x - x3) / (x4 - x3)).powi(2)),
            );
        }
        FTWindow::KaiserBessel => {
            let cen = (x4 + x1) / 2.0;
            let wid = (x4 - x1) / 2.0;
            let arg = (x - cen)
                .mapv(|x| 1.0 - x.powi(2) / wid.powi(2))
                .mapv(|x| x.max(0.0));
            let scale = (bessel_i0::bessel_i0(dx1) - 1.0).max(1e-10);

            fwin = arg.mapv(|x| (bessel_i0::bessel_i0(dx1 * x.sqrt()) - 1.0) / scale);
        }
        FTWindow::Sine => {
            fwin.slice_mut(ndarray::s![i1..=i4]).assign(
                &x.slice(ndarray::s![i1..=i4])
                    .mapv(|x| (std::f64::consts::PI * (x4 - x) / (x4 - x1)).sin()),
            );
        }
        FTWindow::Gaussian => {
            let cen = (x4 + x1) / 2.0;
            fwin = x.mapv(|x| (-(x - cen).powi(2) / (2.0 * dx1.powi(2))).exp());
        }
    }

    Ok(fwin)
}

pub use super::tools::{RebinConfig, RebinMethod, RebinOutput};

/// `(energy, mu, stddev)` arrays returned by [`rebin`].
pub type RebinArrays = (Array1<f64>, Array1<f64>, Array1<f64>);

/// Rebin `(energy, mu)` onto Athena's standard three-region grid.
///
/// Thin ndarray wrapper around [`super::tools::rebin`]; returns
/// `(energy, mu, stddev)`. See [`RebinConfig`] for the region boundaries.
pub fn rebin<T: Into<Array1<f64>>>(
    energy: T,
    mu: T,
    cfg: &RebinConfig,
) -> Result<RebinArrays, super::XAFSError> {
    let energy = DVector::from_vec(energy.into().to_vec());
    let mu = DVector::from_vec(mu.into().to_vec());
    let out = super::tools::rebin(&energy, &mu, cfg)?;
    Ok((
        Array1::from_vec(out.energy.as_slice().to_vec()),
        Array1::from_vec(out.mu.as_slice().to_vec()),
        Array1::from_vec(out.stddev.as_slice().to_vec()),
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::xafs::tests::PARAM_LOADTXT;
    use crate::xafs::tests::TEST_TOL;
    use crate::xafs::tests::TOP_DIR;
    use approx::{assert_abs_diff_eq, assert_abs_diff_ne};
    use data_reader::reader::{load_txt_f64, Delimiter, ReaderParams};
    const ACCEPTABLE_MU_DIFF: f64 = 1e-2;
    const TEST_TOL_FTWINDOW: f64 = 1e-15;

    #[test]
    fn test_smooth() -> Result<(), Box<dyn std::error::Error>> {
        let filepath = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS_smooth.txt";
        let expected_filepath_larch =
            String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS_smooth_larch.txt";
        let xafs_group = io::load_spectrum_QAS_trans(&filepath)?;

        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT)?;
        let expected_data = expected_data.get_col(0);

        let expected_data_larch = load_txt_f64(&expected_filepath_larch, &PARAM_LOADTXT)?;
        let expected_data_larch = expected_data_larch.get_col(0);

        // Convert DVector to Array1 for smooth function
        let x = Array1::from_vec(xafs_group.raw_energy.unwrap().data.as_vec().clone());
        let y = Array1::from_vec(xafs_group.raw_mu.unwrap().data.as_vec().clone());

        let result = smooth(x, y, None, None, None, None, ConvolveForm::Lorentzian)?;

        result
            .iter()
            .zip(expected_data)
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL));

        result
            .iter()
            .zip(expected_data_larch.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = ACCEPTABLE_MU_DIFF));

        Ok(())
    }

    #[test]
    fn test_remove_dups() {
        let arr = DVector::from_vec(vec![0.0, 1.1, 2.2, 2.2, 3.3]);
        let arr = remove_dups(&arr, None, None, None);
        let expected = DVector::from_vec(vec![0.0, 1.1, 2.2, 2.2000001, 3.3]);

        arr.iter().zip(expected.iter()).for_each(|(a, b)| {
            assert_abs_diff_eq!(a, b, epsilon = TEST_TOL);
        });
    }

    #[test]
    fn test_remove_dups_sort() {
        let arr = DVector::from_vec(vec![0.0, 1.1, 2.2, 3.3, 2.2]);
        let arr = remove_dups(&arr, None, None, Some(true));
        let expected = DVector::from_vec(vec![0.0, 1.1, 2.2, 2.2000001, 3.3]);

        arr.iter().zip(expected.iter()).for_each(|(a, b)| {
            assert_abs_diff_eq!(a, b, epsilon = TEST_TOL);
        });
    }

    #[test]
    fn test_remove_dups_unsorted() {
        let arr = DVector::from_vec(vec![0.0, 1.1, 2.2, 3.3, 2.2]);
        let arr = remove_dups(&arr, None, None, Some(false));
        let expected = DVector::from_vec(vec![0.0, 1.1, 2.2, 2.2000001, 3.3]);

        assert_ne!(arr, DVector::from_vec(vec![0., 1.1, 2.2, 2.2000001, 3.3]));
    }

    #[test]
    fn test_find_energy_step() {
        let energy = DVector::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0]);
        let step = find_energy_step(&energy, None, None, None);
        assert_eq!(step, 1.0);
    }

    #[test]
    fn test_find_energy_step_neg() {
        let energy = DVector::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0, 2.0]);
        let step = find_energy_step(&energy, None, None, None);
        assert_eq!(step, 0.25);
    }

    #[test]
    fn test_find_energy_step_sort() {
        let energy = DVector::from_vec(vec![0.0, 1.0, 2.0, 3.0, 4.0, 2.0]);
        let step = find_energy_step(&energy, Some(0.), None, Some(true));
        assert_eq!(step, 0.75);
    }

    #[test]
    fn test_find_e0() {
        use crate::xafs::nshare::ToNalgebra;
        let energy: Array1<f64> = Array1::linspace(0.0, 100.0, 1000);
        let mu = &energy.map(|x| (x - 50.0).powi(3) - (x - 50.0).powi(2) + x);
        let energy_dv = energy.into_nalgebra();
        let mu_dv = mu.clone().into_nalgebra();
        let result = find_e0(&energy_dv, &mu_dv);

        assert_abs_diff_eq!(result.unwrap(), 0.4004004004004004, epsilon = TEST_TOL);
    }

    #[allow(non_snake_case)]
    #[test]
    fn test_KTOE() {
        let expected_KTOE = 3.809982110968585;

        assert_abs_diff_eq!(constants::KTOE, expected_KTOE, epsilon = TEST_TOL);
    }

    #[test]
    fn test_ftwindow_hanning() {
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/window_Hanning.txt";
        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT).unwrap();
        let x = expected_data.get_col(0);
        let y_expected = expected_data.get_col(1);

        let y = ftwindow(
            &Array1::from_vec(x),
            None,
            None,
            None,
            None,
            Some(FTWindow::Hanning),
        )
        .unwrap();

        y.iter()
            .zip(y_expected.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL));
    }
    #[test]
    fn test_ftwindow_parzen() {
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/window_Parzen.txt";
        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT).unwrap();
        let x = expected_data.get_col(0);
        let y_expected = expected_data.get_col(1);

        let y = ftwindow(
            &Array1::from_vec(x),
            None,
            None,
            None,
            None,
            Some(FTWindow::Parzen),
        )
        .unwrap();

        y.iter()
            .zip(y_expected.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL));
    }
    #[test]
    fn test_ftwindow_welch() {
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/window_Welch.txt";
        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT).unwrap();
        let x = expected_data.get_col(0);
        let y_expected = expected_data.get_col(1);

        let y = ftwindow(
            &Array1::from_vec(x),
            None,
            None,
            None,
            None,
            Some(FTWindow::Welch),
        )
        .unwrap();

        y.iter()
            .zip(y_expected.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL));
    }
    #[test]
    fn test_ftwindow_gaussian() {
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/window_Gaussian.txt";
        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT).unwrap();
        let x = expected_data.get_col(0);
        let y_expected = expected_data.get_col(1);

        let y = ftwindow(
            &Array1::from_vec(x),
            None,
            None,
            None,
            None,
            Some(FTWindow::Gaussian),
        )
        .unwrap();

        y.iter()
            .zip(y_expected.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL));
    }
    #[test]
    fn test_ftwindow_sine() {
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/window_Sine.txt";
        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT).unwrap();
        let x = expected_data.get_col(0);
        let y_expected = expected_data.get_col(1);

        let y = ftwindow(
            &Array1::from_vec(x),
            None,
            None,
            None,
            None,
            Some(FTWindow::Sine),
        )
        .unwrap();

        y.iter()
            .zip(y_expected.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL));
    }

    #[test]
    fn test_ftwindow_kaiserbessel() {
        let expected_filepath = String::from(TOP_DIR) + "/tests/testfiles/window_Kaiser-Bessel.txt";
        let expected_data = load_txt_f64(&expected_filepath, &PARAM_LOADTXT).unwrap();
        let x = expected_data.get_col(0);
        let y_expected = expected_data.get_col(1);

        let y = ftwindow(
            &Array1::from_vec(x),
            None,
            None,
            None,
            None,
            Some(FTWindow::KaiserBessel),
        )
        .unwrap();

        y.iter()
            .zip(y_expected.iter())
            .for_each(|(a, b)| assert_abs_diff_eq!(a, &b, epsilon = TEST_TOL_FTWINDOW));
    }
}