rexafs 0.1.2

#![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
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;

/// Trait for xafs utilities
/// functions for f64, Vec<f64>, and ArrayBase<OwnedRepr<f64>, Ix1>
pub trait XAFSUtils {
    fn etok(&self) -> Self;
    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)
    }
}

#[derive(Debug, Clone, Copy, Default)]
pub enum ConvolveForm {
    #[default]
    Lorentzian,
    Gaussian,
    Voigt,
}
/// Smooth a funtion y(x) by convoluting with a lorentzian, gaussian, or voigt function.
///
/// The function is sampled at intervals xstep, and the convolution is performed
/// using FFT convolution. The function is padded with npad points on each side
/// before convolution. The function is interpolated onto a uniform grid with
/// spacing xstep before convolution. The function is returned on the original
/// grid.
///
/// # Arguments
/// * `x` - x values of the function
/// * `y` - y values of the function
/// * `sigma` - primary width parameter for convolving function (default: 1.0)
/// * `gamma` - secondary width parameter for convolving function (default: sigma)
/// * `xstep` - step size for uniform grid onto which the function is interpolated (default: min(x.di/ff())
/// * `npad` - number of points to pad onto each side of the function before convolution (default: 5)/
/// * `conv_form` - form of the convolving function (default: lorentzian)
///
/// # Returns
/// * Result<Array1<f64>, Box<dyn Error>> - smoothed function
///
/// # Example
/// ```
/// use ndarray::Array1;
/// use rexafs::xafs::xafsutils::{smooth, ConvolveForm};
///
/// let x: Array1<f64> = Array1::range(0.0, 10.0, 1.0);
/// let y: Array1<f64> = Array1::range(0.0, 10.0, 1.0);
///
/// let result = smooth(x, y, None, None, None, None, ConvolveForm::Lorentzian);
/// ```
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)
}

/// Function to remove duplicated successive values of an array that is expected to be monotonically increasing.
///
/// For repeated value, the second encountered occurrence (at index i) will be increased by an amount that is the larget of:
/// 1. tiny (default 1e-7)
/// 2. frac (default 1e-6) times the difference between the previous and next values.
///
/// # Arguments
/// * `arr` - Array of values to be checked for duplicates
/// * `tiny` - Minimum value to be added to a duplicate value (default 1e-7)
/// * `frac` - Fraction of the difference between the previous and next values to be added to a duplicate value (default 1e-6)
///
/// # Returns
/// * `arr` - Array with duplicates removed
///
/// # Example
/// ```
/// use rexafs::xafs::xafsutils::remove_dups_array1;
/// use ndarray::Array1;
///
/// let arr = Array1::from_vec(vec![0.0, 1.1, 2.2, 2.2, 3.3]);
/// let arr = remove_dups_array1(arr, None, None, None);
/// assert_eq!(arr, Array1::from_vec(vec![0., 1.1, 2.2, 2.2000001, 3.3]));
/// ```
#[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
}

/// Function to remove duplicated successive values of a DVector that is expected to be monotonically increasing.
///
/// For repeated value, the second encountered occurrence (at index i) will be increased by an amount that is the largest of:
/// 1. tiny (default 1e-7)
/// 2. frac (default 1e-6) times the difference between the previous and next values.
///
/// # Arguments
/// * `arr` - DVector of values to be checked for duplicates
/// * `tiny` - Minimum value to be added to a duplicate value (default 1e-7)
/// * `frac` - Fraction of the difference between the previous and next values to be added to a duplicate value (default 1e-6)
/// * `sort` - Sort the vector before removing duplicates (default false)
///
/// # Returns
/// * `DVector<f64>` - Vector with duplicates removed
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
}

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())
}

/// Function to find the energy step of an array of energies.
/// It ignores the smallest fraction of energy steps (frac_ignore) and then averages the next nave steps.
///
/// # Arguments
/// * `energy` - Array of energies
/// * `frac_ignore` - Fraction of energy steps to ignore (default 0.01)
/// * `nave` - Number of energy steps to average (default 10)
/// * `sort` - Sort the array before finding the energy step (default false)
///
/// # Returns
/// * `estep` - Average energy step
///
/// # Example
/// ```
/// use rexafs::xafs::xafsutils::find_energy_step_array1;
/// use ndarray::array;
///
/// let energy = array![0.0, 1.1, 2.2, 2.2, 3.3];
/// let estep = find_energy_step_array1(energy, None, None, None);
/// assert_eq!(estep, 0.7333333333333333);
/// ```
#[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;
}
/// Calculate the $E_0$, the energy threshold of absoption, or the edge energy, given $\mu(E)$.
///
/// $E_0$ is found as the point with maximum derivative with some checks to avoid spurious glitches.
///
/// # Arguments
/// * `energy` - Array of energies
/// * `mu` - Array of absorption coefficients
///
/// # Returns
/// Result<e0: f64, Box<dyn Error>>
/// * `e0` - Energy threshold of absoption, or the edge energy
///
/// # Example
/// ```
/// use rexafs::xafs::xafsutils::find_e0_array1;
/// use ndarray::Array1;
///
/// 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 result = find_e0_array1(energy.clone(), mu.clone());
/// assert_eq!(result.unwrap(), 0.4004004004004004);
///
/// // Result calculated by Larch is 0.3003003003003003
/// ```

#[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)
}

/// Internal function used for find_e0.
///
/// # Arguments
/// * `energy` - Array of energies
/// * `mu` - Array of absorption coefficients
/// * `estep` - Energy step (default: find_energy_step(energy)/2.0)
/// * `use_smooth` - Use smoothed derivative (default: false)
///
/// # Returns
/// Result<(e0: f64, imax: usize, estep: f64), Box<dyn Error>>
/// * `e0` - Energy threshold of absoption, or the edge energy
/// * `imax` - Index of maximum derivative
/// * `estep` - Energy step
///
/// # Example
/// ```
/// use rexafs::xafs::xafsutils::_find_e0_array1;
/// use ndarray::Array1;
///
/// 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 result = _find_e0_array1(energy.clone(), mu.clone(), None, None);
/// assert_eq!(result.unwrap(), (1.001001001001001, 10, 0.05005005005004648));
///
/// // the result obtained by xraylarch is (1.001001001001001, 10, 0.05005005005004648)
/// ```
#[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))
}

/// DVector version: Find the energy step of an array of energies.
/// It ignores the smallest fraction of energy steps (frac_ignore) and then averages the next nave steps.
///
/// # Arguments
/// * `energy` - DVector of energies
/// * `frac_ignore` - Fraction of energy steps to ignore (default 0.01)
/// * `nave` - Number of energy steps to average (default 10)
/// * `sort` - Sort the array before finding the energy step (default false)
///
/// # Returns
/// * `estep` - Average energy step
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
}

/// DVector version: Internal function used for find_e0.
///
/// # Arguments
/// * `energy` - DVector of energies
/// * `mu` - DVector of absorption coefficients
/// * `estep` - Energy step (default: find_energy_step(energy)/2.0)
/// * `use_smooth` - Use smoothed derivative (default: false)
///
/// # Returns
/// Result<(e0: f64, imax: usize, estep: f64), Box<dyn Error>>
/// * `e0` - Energy threshold of absorption, or the edge energy
/// * `imax` - Index of maximum derivative
/// * `estep` - Energy step
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))
}

/// DVector version: Calculate the E_0, the energy threshold of absorption, or the edge energy, given μ(E).
///
/// E_0 is found as the point with maximum derivative with some checks to avoid spurious glitches.
///
/// # Arguments
/// * `energy` - DVector of energies
/// * `mu` - DVector of absorption coefficients
///
/// # Returns
/// Result<e0: f64, Box<dyn Error>>
/// * `e0` - Energy threshold of absorption, or the edge energy
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)
}

#[derive(Debug, Clone, Copy, Default, PartialEq, Serialize, Deserialize)]
pub enum FTWindow {
    #[default]
    Hanning, // Hanning window, cosine-squared tamper
    Parzen,       // Parzen window, linear tamper
    Welch,        // Welch window, quadratic tamper
    Gaussian,     // Gaussian window, Gaussian (normal) tamper
    Sine,         // Sine window, sine function window
    KaiserBessel, // Kaiser-Bessel function-derived window
    FHanning,     // I am not sure what this is. It is in the Larch code, but it is not used.
}

impl FTWindow {
    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))
    }
}

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));
    }
}