rexafs 0.1.0

// Standard library dependencies

// External dependencies
use easyfft::prelude::{DynRealFft, DynRealIfft};
use easyfft::{dyn_size::realfft::DynRealDft, num_complex::Complex};
use nalgebra::{DVector, Owned};
use ndarray::{
    Array, Array1, ArrayBase, ArrayView, ArrayView1, Axis, Ix, Ix1, OwnedRepr, ViewRepr,
};
use num_complex::Complex64;
use serde::{Deserialize, Serialize};

use derivative::Derivative;

// load dependencies
use super::errors::FFTError;

// Load local traits
use super::mathutils::MathUtils;
use super::xafsutils::ftwindow;
use crate::xafs::xafsutils::FTWindow;

#[derive(Derivative, Debug, Clone, Serialize, Deserialize)]
#[derivative(PartialEq)]
#[serde(default)]
pub struct XrayFFTF {
    pub rmax_out: Option<f64>,
    pub window: Option<FTWindow>,
    pub dk: Option<f64>,
    pub dk2: Option<f64>,
    pub kmin: Option<f64>,
    pub kmax: Option<f64>,
    pub kweight: Option<f64>,
    pub nfft: Option<usize>,
    pub kstep: Option<f64>,
    pub r: Option<ArrayBase<OwnedRepr<f64>, Ix1>>,
    // currently asking for serde support in the easyfft crate
    #[derivative(PartialEq = "ignore")]
    pub chir: Option<DynRealDft<f64>>,
    pub chir_mag: Option<ArrayBase<OwnedRepr<f64>, Ix1>>,
    pub kwin: Option<ArrayBase<OwnedRepr<f64>, Ix1>>,
}

impl Default for XrayFFTF {
    fn default() -> Self {
        XrayFFTF {
            rmax_out: Some(10.0),
            window: Some(FTWindow::KaiserBessel),
            dk: Some(1.),
            dk2: None,
            kmin: Some(2.0),
            kmax: Some(15.0),
            kweight: Some(2.0),
            nfft: Some(2048),
            kstep: None,
            r: None,
            chir: None,
            chir_mag: None,
            kwin: None,
        }
    }
}

impl XrayFFTF {
    pub fn new() -> XrayFFTF {
        XrayFFTF::default()
    }

    pub fn fill_parameter(&mut self, k: ArrayBase<ViewRepr<&f64>, Ix1>) -> &mut Self {
        if self.kweight.is_none() {
            self.kweight = Some(2.0);
        }

        self.kweight = Some(self.kweight.unwrap().max(0.0).floor());

        if self.kstep.is_none() {
            self.kstep = Some(k[1] - k[0]);
        }

        if self.kmin.is_none() {
            self.kmin = Some(k[0]);
        }

        if self.kmax.is_none() {
            self.kmax = Some(k[k.len() - 1]);
        }

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

        if self.dk2.is_none() {
            self.dk2 = self.dk;
        }

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

        if self.rmax_out.is_none() {
            self.rmax_out = Some(10.0);
        }

        self
    }

    pub fn xftf_prep(
        &mut self,
        k: ArrayBase<ViewRepr<&f64>, Ix1>,
        chi: ArrayBase<ViewRepr<&f64>, Ix1>,
    ) -> Result<
        (
            ArrayBase<OwnedRepr<f64>, Ix1>,
            ArrayBase<OwnedRepr<f64>, Ix1>,
        ),
        FFTError,
    > {
        if self.nfft.is_some_and(|n| n < 2) {
            return Err(FFTError::InvalidParameter {
                parameter: "nfft".into(),
                reason: "must be at least 2".into(),
            });
        }
        for (name, value, strictly_positive) in [
            ("kstep", self.kstep, true),
            ("kweight", self.kweight, false),
            ("dk", self.dk, false),
            ("dk2", self.dk2, false),
            ("rmax_out", self.rmax_out, false),
        ] {
            if value.is_some_and(|v| !v.is_finite() || v < 0.0 || (strictly_positive && v == 0.0)) {
                return Err(FFTError::InvalidParameter {
                    parameter: name.into(),
                    reason: "must be finite and nonnegative (kstep must be positive)".into(),
                });
            }
        }
        if self.kmin.is_some_and(|v| !v.is_finite())
            || self.kmax.is_some_and(|v| !v.is_finite())
            || self.kmin.zip(self.kmax).is_some_and(|(lo, hi)| lo >= hi)
        {
            return Err(FFTError::InvalidParameter {
                parameter: "kmin/kmax".into(),
                reason: "must be finite with kmin < kmax".into(),
            });
        }
        if k.len() != chi.len() {
            return Err(FFTError::InterpolationFailed {
                reason: "k/chi length mismatch".to_string(),
            });
        }
        if k.len() < 2 {
            let (kmin, kmax) = if k.is_empty() {
                (0.0, 0.0)
            } else {
                (k[0], k[k.len() - 1])
            };
            return Err(FFTError::InsufficientPoints {
                min: 2,
                actual: k.len(),
                kmin,
                kmax,
            });
        }

        self.fill_parameter(k);
        let kweight = self.kweight.unwrap() as i32;
        let k_max = k.iter().max_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
        let npts = (1.01 + k_max / self.kstep.unwrap()) as usize;
        let k_max = k_max.max(self.kmax.unwrap() + self.dk2.unwrap());
        let k_ = Array1::range(0.0, k_max + self.kstep.unwrap(), self.kstep.unwrap());

        let chi_ = if let (Some(k_slice), Some(chi_slice)) = (k.as_slice(), chi.as_slice()) {
            k_.interpolate(k_slice, chi_slice)
        } else {
            k_.interpolate(&k.to_vec(), &chi.to_vec())
        }
        .map_err(|e| FFTError::InterpolationFailed {
            reason: e.to_string(),
        })?;
        let win = self
            .window
            .unwrap()
            .window(&k_, self.kmin, self.kmax, self.dk, self.dk2)
            .map_err(|e| FFTError::WindowCalculationFailed {
                reason: e.to_string(),
            })?;
        let win = (win).slice_axis(Axis(0), (0..npts).into()).to_owned();
        let chi_ = &chi_.slice_axis(Axis(0), (0..npts).into())
            * &k_
                .slice_axis(Axis(0), (0..npts).into())
                .map(|x| x.powi(kweight));

        Ok((chi_, win))
    }

    pub fn xftf(
        &mut self,
        k: ArrayBase<ViewRepr<&f64>, Ix1>,
        chi: ArrayBase<ViewRepr<&f64>, Ix1>,
    ) -> Result<&mut Self, FFTError> {
        let (cchi, win) = self.xftf_prep(k, chi)?;

        let cchi_fft = xftf_fast(cchi.view(), self.nfft.unwrap(), self.kstep.unwrap());

        let rstep = std::f64::consts::PI / self.kstep.unwrap() / self.nfft.unwrap() as f64;

        // The length of r is different by 1 between xraylarch and rexafs. This is due to the implementation of FFT.
        let irmax =
            (self.nfft.unwrap() / 2 + 1).min((1.01 + self.rmax_out.unwrap() / rstep) as usize);

        self.r = Some(Array1::range(0.0, irmax as f64 * rstep, rstep));

        self.chir = Some(cchi_fft.clone());
        self.chir_mag = Some(cchi_fft[0..irmax].norm());
        self.kwin = Some(win);

        Ok(self)
    }

    pub fn get_rmax_out(&self) -> Option<&f64> {
        self.rmax_out.as_ref()
    }

    pub fn get_window(&self) -> Option<&FTWindow> {
        self.window.as_ref()
    }

    pub fn get_dk(&self) -> Option<&f64> {
        self.dk.as_ref()
    }

    pub fn get_dk2(&self) -> Option<&f64> {
        self.dk2.as_ref()
    }

    pub fn get_kmin(&self) -> Option<&f64> {
        self.kmin.as_ref()
    }

    pub fn get_kmax(&self) -> Option<&f64> {
        self.kmax.as_ref()
    }

    pub fn get_kweight(&self) -> Option<&f64> {
        self.kweight.as_ref()
    }

    pub fn get_r(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        Some(self.r.as_ref()?.view())
    }

    pub fn get_chir(&self) -> Option<&DynRealDft<f64>> {
        self.chir.as_ref()
    }

    pub fn get_chir_real(&self) -> Option<ArrayBase<OwnedRepr<f64>, Ix1>> {
        let len_r = self.r.as_ref()?.len();
        let chir = self.chir.as_ref()?;
        Some(Array1::from_iter(chir.iter().take(len_r).map(|x| x.re)))
    }

    pub fn get_chir_imag(&self) -> Option<ArrayBase<OwnedRepr<f64>, Ix1>> {
        let len_r = self.r.as_ref()?.len();
        let chir = self.chir.as_ref()?;
        Some(Array1::from_iter(chir.iter().take(len_r).map(|x| x.im)))
    }

    pub fn get_chir_mag(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        Some(self.chir_mag.as_ref()?.view())
    }
    pub fn get_kwin(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        Some(self.kwin.as_ref()?.view())
    }

    pub fn get_kstep(&self) -> Option<&f64> {
        self.kstep.as_ref()
    }
}

#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
#[serde(default)]
pub struct XrayFFTR {
    pub qmax_out: Option<f64>,
    pub window: Option<FTWindow>,
    pub dr: Option<f64>,
    pub dr2: Option<f64>,
    pub rmin: Option<f64>,
    pub rmax: Option<f64>,
    pub rweight: Option<f64>,
    pub nfft: Option<usize>,
    pub kstep: Option<f64>,
    pub q: Option<ArrayBase<OwnedRepr<f64>, Ix1>>,
    pub chiq: Option<ArrayBase<OwnedRepr<f64>, Ix1>>,
    pub rwin: Option<ArrayBase<OwnedRepr<f64>, Ix1>>,
}

impl Default for XrayFFTR {
    fn default() -> Self {
        XrayFFTR {
            qmax_out: Some(10.0),
            window: Some(FTWindow::KaiserBessel),
            dr: Some(1.),
            dr2: None,
            rmin: Some(0.),
            rmax: Some(20.0),
            rweight: Some(0.0),
            nfft: Some(2048),
            kstep: None,
            q: None,
            chiq: None,
            rwin: None,
        }
    }
}

impl XrayFFTR {
    pub fn new() -> XrayFFTR {
        XrayFFTR::default()
    }

    pub fn fill_parameter(&mut self, r: ArrayBase<ViewRepr<&f64>, Ix1>) -> &mut Self {
        if self.rweight.is_none() {
            self.rweight = Some(0.0);
        }

        self.rweight = Some(self.rweight.unwrap().max(0.0).floor());

        if self.rmin.is_none() {
            self.rmin = Some(r[0]);
        }

        if self.rmax.is_none() {
            self.rmax = Some(r[r.len() - 1]);
        }

        if self.dr.is_none() {
            self.dr = Some(1.0);
        }

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

        if self.qmax_out.is_none() {
            self.qmax_out = Some(10.0);
        }

        if self.kstep.is_none() {
            self.kstep = Some(std::f64::consts::PI / (r[1] - r[0]) / self.nfft.unwrap() as f64);
        }
        self
    }

    pub fn xftr_prep(
        &mut self,
        r: ArrayBase<ViewRepr<&f64>, Ix1>,
        chir: &DynRealDft<f64>,
    ) -> Result<(DynRealDft<f64>, ArrayBase<OwnedRepr<f64>, Ix1>), FFTError> {
        self.fill_parameter(r);
        let rweight = self.rweight.unwrap() as i32;
        let nfft = self.nfft.unwrap();
        let r_len = chir.len();
        let rstep = std::f64::consts::PI / self.kstep.unwrap() / nfft as f64;

        let r_ = Array1::range(0.0, r_len as f64 * rstep, rstep);

        let win = if rweight == 0 {
            ftwindow(&r_, self.rmin, self.rmax, self.dr, self.dr2, self.window).map_err(|e| {
                FFTError::WindowCalculationFailed {
                    reason: e.to_string(),
                }
            })?
        } else {
            ftwindow(&r_, self.rmin, self.rmax, self.dr, self.dr2, self.window).map_err(|e| {
                FFTError::WindowCalculationFailed {
                    reason: e.to_string(),
                }
            })? * &r_.map(|x| x.powi(rweight))
        };

        let chir_win = chir
            .iter()
            .zip(win.iter())
            .map(|(x, y)| x * y)
            .collect::<Vec<Complex<f64>>>();

        let chir_win = DynRealDft::new(*chir.get_offset(), &chir_win[1..], nfft);

        Ok((chir_win, win))
    }

    pub fn xftr(
        &mut self,
        r: ArrayBase<ViewRepr<&f64>, Ix1>,
        chir: &DynRealDft<f64>,
    ) -> Result<&mut Self, FFTError> {
        let (chir_win, win) = self.xftr_prep(r, chir)?;
        let nfft = self.nfft.unwrap();
        let out = xftr_fast(&chir_win, nfft, self.kstep.unwrap());

        let q = Array1::linspace(
            0.0,
            self.qmax_out.unwrap(),
            (1.05 + self.qmax_out.unwrap() / self.kstep.unwrap()) as usize,
        );

        self.q = Some(q);
        self.rwin = Some(win);
        self.chiq = Some(out);

        Ok(self)
    }

    pub fn get_q(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        Some(self.q.as_ref()?.view())
    }

    pub fn get_chiq(&self) -> Option<ArrayBase<OwnedRepr<f64>, Ix1>> {
        let len_q = self.q.as_ref()?.len();

        Some(
            self.chiq
                .clone()?
                .slice_axis(Axis(0), (0..len_q).into())
                .to_owned(),
        )
    }

    pub fn get_rwin(&self) -> Option<ArrayBase<ViewRepr<&f64>, Ix1>> {
        Some(self.rwin.as_ref()?.view())
    }

    pub fn get_kstep(&self) -> Option<&f64> {
        self.kstep.as_ref()
    }

    pub fn get_rweight(&self) -> Option<&f64> {
        self.rweight.as_ref()
    }

    pub fn get_nfft(&self) -> Option<&usize> {
        self.nfft.as_ref()
    }

    pub fn get_window(&self) -> Option<&FTWindow> {
        self.window.as_ref()
    }
}

pub fn xftf_fast(chi: ArrayBase<ViewRepr<&f64>, Ix1>, nfft: usize, kstep: f64) -> DynRealDft<f64> {
    let mut cchi = vec![0.0_f64; nfft];
    for (output, input) in cchi.iter_mut().zip(chi.iter()) {
        *output = *input;
    }

    let mut freq = cchi.real_fft();

    freq *= kstep / (std::f64::consts::PI).sqrt();

    freq
}

pub fn xftr_fast(
    chir: &DynRealDft<f64>,
    nfft: usize,
    kstep: f64,
) -> ArrayBase<OwnedRepr<f64>, Ix1> {
    let cchi = if chir.len() < nfft / 2 + 1 {
        let mut freq_bin = vec![Complex::new(0.0, 0.0); nfft - 1];
        freq_bin[..chir.len() - 1].copy_from_slice(chir.get_frequency_bins());
        DynRealDft::new(*chir.get_offset(), &freq_bin, nfft)
    } else {
        chir.clone()
    };

    let mut chi = Array1::from(cchi.real_ifft());

    chi *= std::f64::consts::PI.sqrt() / kstep / nfft as f64;

    chi
}

pub fn xftf_fast_nalgebra(chi: &DVector<f64>, nfft: usize, kstep: f64) -> DynRealDft<f64> {
    let mut cchi = vec![0.0_f64; nfft];
    let length = chi.len().min(nfft);
    cchi[..length].copy_from_slice(&chi.as_slice()[..length]);

    let mut freq = cchi.real_fft();

    freq *= kstep / std::f64::consts::PI.sqrt();

    freq
}

pub fn xftr_fast_nalgebra(chir: &DynRealDft<f64>, nfft: usize, kstep: f64) -> DVector<f64> {
    let cchi = if chir.len() < nfft / 2 + 1 {
        let mut freq_bin = vec![Complex::new(0.0, 0.0); nfft - 1];
        freq_bin[..chir.len() - 1].copy_from_slice(chir.get_frequency_bins());
        DynRealDft::new(*chir.get_offset(), &freq_bin, nfft)
    } else {
        chir.clone()
    };

    let mut chi = DVector::from(cchi.real_ifft().to_vec());

    chi *= (std::f64::consts::PI).sqrt() / kstep / nfft as f64;

    chi
}

pub trait XFFT {
    fn xftf_fast(&self, nfft: usize, kstep: f64) -> DynRealDft<f64>;
}

impl XFFT for ArrayBase<OwnedRepr<f64>, Ix1> {
    fn xftf_fast(&self, nfft: usize, kstep: f64) -> DynRealDft<f64> {
        xftf_fast(self.view(), nfft, kstep)
    }
}

impl XFFT for ArrayBase<ViewRepr<&f64>, Ix1> {
    fn xftf_fast(&self, nfft: usize, kstep: f64) -> DynRealDft<f64> {
        xftf_fast(*self, nfft, kstep)
    }
}

impl XFFT for DVector<f64> {
    fn xftf_fast(&self, nfft: usize, kstep: f64) -> DynRealDft<f64> {
        xftf_fast_nalgebra(self, nfft, kstep)
    }
}

pub trait XFFTReverse<T> {
    fn xftr_fast(&self, nfft: usize, kstep: f64) -> T;
}

impl XFFTReverse<ArrayBase<OwnedRepr<f64>, Ix1>> for DynRealDft<f64> {
    fn xftr_fast(&self, nfft: usize, kstep: f64) -> ArrayBase<OwnedRepr<f64>, Ix1> {
        xftr_fast(self, nfft, kstep)
    }
}

impl XFFTReverse<DVector<f64>> for DynRealDft<f64> {
    fn xftr_fast(&self, nfft: usize, kstep: f64) -> DVector<f64> {
        xftr_fast_nalgebra(self, nfft, kstep)
    }
}

pub trait FFTUtils<T> {
    fn realimg(&self) -> T;
    fn re(&self) -> T;
    fn im(&self) -> T;
    fn norm(&self) -> T;
    fn norm_sqr(&self) -> T;
}

impl FFTUtils<ArrayBase<OwnedRepr<f64>, Ix1>> for DynRealDft<f64> {
    fn realimg(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().flat_map(|x| vec![x.re, x.im]))
    }

    fn re(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.re))
    }

    fn im(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.im))
    }

    fn norm(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.norm()))
    }

    fn norm_sqr(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.norm_sqr()))
    }
}

impl FFTUtils<DVector<f64>> for DynRealDft<f64> {
    fn realimg(&self) -> DVector<f64> {
        DVector::from_iterator(self.len() * 2, self.iter().flat_map(|x| vec![x.re, x.im]))
    }

    fn re(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.re))
    }

    fn im(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.im))
    }

    fn norm(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.norm()))
    }

    fn norm_sqr(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.norm_sqr()))
    }
}

impl FFTUtils<ArrayBase<OwnedRepr<f64>, Ix1>> for [Complex<f64>] {
    fn realimg(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().flat_map(|x| vec![x.re, x.im]))
    }

    fn re(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.re))
    }

    fn im(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.im))
    }

    fn norm(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.norm()))
    }

    fn norm_sqr(&self) -> Array1<f64> {
        Array1::from_iter(self.iter().map(|x| x.norm_sqr()))
    }
}

impl FFTUtils<DVector<f64>> for [Complex<f64>] {
    fn realimg(&self) -> DVector<f64> {
        DVector::from_iterator(self.len() * 2, self.iter().flat_map(|x| vec![x.re, x.im]))
    }

    fn re(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.re))
    }

    fn im(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.im))
    }

    fn norm(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.norm()))
    }

    fn norm_sqr(&self) -> DVector<f64> {
        DVector::from_iterator(self.len(), self.iter().map(|x| x.norm_sqr()))
    }
}

// #[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
// pub struct DynRealDftWrapper {
//     zeroth_bin: f64,

//     frequency_bins: Vec<Complex<f64>>,

//     original_length: usize,
// }

// impl From<DynRealDft<f64>> for DynRealDftWrapper {
//     fn from(dft: DynRealDft<f64>) -> Self {
//         DynRealDftWrapper {
//             zeroth_bin: dft.get_offset().clone(),
//             frequency_bins: dft.get_frequency_bins().to_vec(),
//             original_length: dft.get_original_length(),
//         }
//     }
// }

// impl PartialEq for DynRealDft<f64> {
//     fn eq(&self, other: &Self) -> bool {
//         self.len() == other.len()
//             && self.get_offset() == other.get_offset()
//             && self.get_frequency_bins() == other.get_frequency_bins()
//     }
// }

// #[derive(Debug, Clone, PartialEq, Default)]
// pub struct XrayFFTR {}

#[cfg(test)]
mod tests {
    use easyfft::prelude::*;
    use ndarray::Array1;

    use super::*;
    use crate::xafs::io;
    use crate::xafs::nshare::ToNalgebra;
    use approx::{assert_abs_diff_eq, assert_relative_eq};

    use data_reader::reader::{load_txt_f64, Delimiter, ReaderParams};

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

    use crate::xafs::background::BackgroundMethod;
    use crate::xafs::background::AUTOBK;

    const ACCEPTABLE_MU_DIFF: f64 = 1e-6;
    const CHI_MSE_TOL: f64 = 1e-2;
    const CHI_Q_TOL: f64 = 1e-1;

    #[test]
    fn test_xftf_fast() {
        let x: Array1<f64> = Array1::linspace(0., 10., 10);
        let sin_x = x.map(|x| x.sin());
        let nfft = 10;
        let kstep = 1.;
        let fft = xftf_fast(sin_x.view(), nfft, kstep);

        let norm: DVector<f64> = fft.norm();

        let expected_norm = DVector::from(vec![
            0.6822515304148188,
            0.999816632055004,
            2.4133321684349966,
            0.35447122637608214,
            0.16620199767343982,
            0.1252841340812192,
        ]);

        assert_abs_diff_eq!(norm, expected_norm, epsilon = 1e-16);
    }

    #[test]
    fn test_xftr_fast() {
        let x: Array1<f64> = Array1::linspace(0., 10., 1024);
        let sin_x = x.map(|x| x.sin());
        let nfft = 1024;
        let kstep = 0.1;

        let fft = xftf_fast(sin_x.view(), nfft, kstep);
        let ifft = xftr_fast(&fft, nfft, kstep);

        sin_x.iter().zip(ifft.iter()).for_each(|(x, y)| {
            assert_abs_diff_eq!(x, y, epsilon = 1e-12);
        });
    }

    #[test]
    fn test_xftf_fast_nalgebra() {
        let x: DVector<f64> = Array1::linspace(0., 10., 10).into_nalgebra();
        let sin_x = x.map(|x| x.sin());
        let nfft = 10;
        let kstep = 1.;
        let fft = xftf_fast_nalgebra(&sin_x, nfft, kstep);

        let norm: DVector<f64> = fft.norm();

        let expected_norm = DVector::from(vec![
            0.6822515304148188,
            0.999816632055004,
            2.4133321684349966,
            0.35447122637608214,
            0.16620199767343982,
            0.1252841340812192,
        ]);

        assert_abs_diff_eq!(norm, expected_norm, epsilon = 1e-16);
    }

    #[test]
    fn test_xftr_fast_nalgebra() {
        let x: DVector<f64> = Array1::linspace(0., 10., 1024).into_nalgebra();
        let sin_x = x.map(|x| x.sin());
        let nfft = 1024;
        let kstep = 0.1;

        let fft = xftf_fast_nalgebra(&sin_x, nfft, kstep);
        let ifft = xftr_fast_nalgebra(&fft, nfft, kstep);

        sin_x.iter().zip(ifft.iter()).for_each(|(x, y)| {
            assert_abs_diff_eq!(x, y, epsilon = 1e-12);
        });
    }

    #[test]
    #[allow(non_snake_case)]
    fn test_Xray_FFTF() -> Result<(), Box<dyn std::error::Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut xafs_test_group = io::load_spectrum_QAS_trans(&path).unwrap();

        xafs_test_group.set_background_method(Some(BackgroundMethod::AUTOBK(AUTOBK {
            rbkg: Some(1.4),
            kweight: Some(2),
            ..Default::default()
        })))?;
        xafs_test_group.calc_background()?;

        xafs_test_group.xftf = Some(XrayFFTF {
            window: Some(FTWindow::Hanning),
            dk: Some(1.0),
            kmin: Some(2.0),
            kmax: Some(15.0),
            kweight: Some(2.0),
            ..Default::default()
        });
        xafs_test_group.fft()?;

        let larch_r_path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS_xftf_larch.txt";
        let larch_r = load_txt_f64(&larch_r_path, &PARAM_LOADTXT).unwrap();

        let r_expected = larch_r.get_col(0);
        let chir_expected = larch_r.get_col(1);

        let r = xafs_test_group.r().unwrap();
        let chir = xafs_test_group.chir_mag().unwrap();

        r.iter().zip(r_expected.iter()).for_each(|(x, y)| {
            assert_abs_diff_eq!(x, y, epsilon = TEST_TOL);
        });

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

        // println!("mse: {}", mse);
        assert!(mse < CHI_MSE_TOL);

        // chir.iter().zip(chir_expected.iter()).for_each(|(x, y)| {
        //     println!("[{}, {}],", x, y);
        // });

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

        // println!("mse: {}", mse);

        // assert!(mse < CHI_MSE_TOL);

        Ok(())
    }

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

        xafs_test_group.set_background_method(Some(BackgroundMethod::AUTOBK(AUTOBK {
            rbkg: Some(1.4),
            kweight: Some(2),
            ..Default::default()
        })))?;
        xafs_test_group.calc_background()?;

        xafs_test_group.xftf = Some(XrayFFTF {
            window: Some(FTWindow::Hanning),
            dk: Some(1.0),
            kmin: Some(2.0),
            kmax: Some(15.0),
            kweight: Some(2.0),
            ..Default::default()
        });
        xafs_test_group.fft()?;

        let chir_real = xafs_test_group.chir_real().unwrap();
        let chir_imag = xafs_test_group.chir_imag().unwrap();
        assert_eq!(chir_real.len(), chir_imag.len());

        let has_distinct_value = chir_real
            .iter()
            .zip(chir_imag.iter())
            .any(|(re, im)| (re - im).abs() > TEST_TOL);
        assert!(has_distinct_value);

        Ok(())
    }

    #[test]
    #[allow(non_snake_case)]
    fn test_XrayFFTR() -> Result<(), Box<dyn std::error::Error>> {
        let path = String::from(TOP_DIR) + "/tests/testfiles/Ru_QAS.dat";
        let mut xafs_test_group = io::load_spectrum_QAS_trans(&path).unwrap();

        xafs_test_group.set_background_method(Some(BackgroundMethod::AUTOBK(AUTOBK {
            rbkg: Some(1.4),
            ..Default::default()
        })))?;
        xafs_test_group.calc_background()?;

        xafs_test_group.xftf = Some(XrayFFTF {
            window: Some(FTWindow::Hanning),
            dk: Some(f64::EPSILON),
            kmin: Some(0.0),
            kmax: Some(15.0),
            kweight: Some(2.0),
            ..Default::default()
        });
        xafs_test_group.fft()?;

        xafs_test_group.xftr = Some(XrayFFTR {
            window: Some(FTWindow::Hanning),
            rweight: Some(0.0),
            dr: Some(f64::EPSILON),
            rmin: Some(0.0),
            rmax: Some(10.0),
            ..Default::default()
        });
        xafs_test_group.ifft()?;

        let q = xafs_test_group.q().unwrap();
        let chiq = xafs_test_group.chiq().unwrap();

        // println!("q: {:?}", q.len());
        // println!("chiq: {:?}", chiq.len());

        assert!(q.len() == chiq.len());

        // Convert DVector to Array1 for multiplication with Array1 view
        let chi_dv = xafs_test_group.chi_kweighted().unwrap();
        let chi_array = Array1::from_vec(chi_dv.data.as_vec().clone());
        let chi = chi_array * xafs_test_group.xftf.unwrap().get_kwin().unwrap();

        let chi = chi.slice_axis(Axis(0), (0..chi.len()).into()).to_owned();

        chi.iter().zip(chiq.iter()).for_each(|(x, y)| {
            assert_relative_eq!(x, y, epsilon = CHI_Q_TOL);
        });

        Ok(())
    }
}