vector_ta/indicators/moving_averages/

dma.rs

#[cfg(all(feature = "python", feature = "cuda"))]
use crate::cuda::cuda_available;
#[cfg(all(feature = "python", feature = "cuda"))]
use crate::cuda::moving_averages::CudaDma;
#[cfg(all(feature = "python", feature = "cuda"))]
use crate::cuda::moving_averages::DeviceArrayF32;
#[cfg(all(feature = "python", feature = "cuda"))]
use cust::context::Context;
#[cfg(all(feature = "python", feature = "cuda"))]
use numpy::PyUntypedArrayMethods;
#[cfg(feature = "python")]
use numpy::{IntoPyArray, PyArray1, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2};
#[cfg(feature = "python")]
use pyo3::exceptions::PyValueError;
#[cfg(feature = "python")]
use pyo3::prelude::*;
#[cfg(feature = "python")]
use pyo3::types::PyDict;
#[cfg(all(feature = "python", feature = "cuda"))]
use std::sync::Arc;

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
use serde::{Deserialize, Serialize};
#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
use wasm_bindgen::prelude::*;

use crate::utilities::data_loader::{source_type, Candles};
use crate::utilities::enums::Kernel;
use crate::utilities::helpers::{
    alloc_with_nan_prefix, detect_best_batch_kernel, detect_best_kernel, init_matrix_prefixes,
    make_uninit_matrix,
};
#[cfg(feature = "python")]
use crate::utilities::kernel_validation::validate_kernel;
use aligned_vec::{AVec, CACHELINE_ALIGN};

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
use core::arch::x86_64::*;

#[cfg(not(target_arch = "wasm32"))]
use rayon::prelude::*;

use std::convert::AsRef;
use std::error::Error;
use std::mem::MaybeUninit;
use thiserror::Error;

impl<'a> AsRef<[f64]> for DmaInput<'a> {
    #[inline(always)]
    fn as_ref(&self) -> &[f64] {
        match &self.data {
            DmaData::Slice(slice) => slice,
            DmaData::Candles { candles, source } => source_type(candles, source),
        }
    }
}

#[derive(Debug, Clone)]
pub enum DmaData<'a> {
    Candles {
        candles: &'a Candles,
        source: &'a str,
    },
    Slice(&'a [f64]),
}

#[derive(Debug, Clone)]
pub struct DmaOutput {
    pub values: Vec<f64>,
}

#[derive(Debug, Clone)]
#[cfg_attr(
    all(target_arch = "wasm32", feature = "wasm"),
    derive(Serialize, Deserialize)
)]
pub struct DmaParams {
    pub hull_length: Option<usize>,
    pub ema_length: Option<usize>,
    pub ema_gain_limit: Option<usize>,
    pub hull_ma_type: Option<String>,
}

impl Default for DmaParams {
    fn default() -> Self {
        Self {
            hull_length: Some(7),
            ema_length: Some(20),
            ema_gain_limit: Some(50),
            hull_ma_type: Some("WMA".to_string()),
        }
    }
}

#[derive(Debug, Clone)]
pub struct DmaInput<'a> {
    pub data: DmaData<'a>,
    pub params: DmaParams,
}

impl<'a> DmaInput<'a> {
    #[inline]
    pub fn from_candles(c: &'a Candles, s: &'a str, p: DmaParams) -> Self {
        Self {
            data: DmaData::Candles {
                candles: c,
                source: s,
            },
            params: p,
        }
    }

    #[inline]
    pub fn from_slice(sl: &'a [f64], p: DmaParams) -> Self {
        Self {
            data: DmaData::Slice(sl),
            params: p,
        }
    }

    #[inline]
    pub fn with_default_candles(c: &'a Candles) -> Self {
        Self::from_candles(c, "close", DmaParams::default())
    }

    #[inline]
    pub fn get_hull_length(&self) -> usize {
        self.params.hull_length.unwrap_or(7)
    }

    #[inline]
    pub fn get_ema_length(&self) -> usize {
        self.params.ema_length.unwrap_or(20)
    }

    #[inline]
    pub fn get_ema_gain_limit(&self) -> usize {
        self.params.ema_gain_limit.unwrap_or(50)
    }

    #[inline]
    pub fn get_hull_ma_type(&self) -> String {
        self.params
            .hull_ma_type
            .clone()
            .unwrap_or_else(|| "WMA".to_string())
    }

    #[inline]
    pub fn hull_ma_type_str(&self) -> &str {
        self.params.hull_ma_type.as_deref().unwrap_or("WMA")
    }
}

#[derive(Clone, Debug)]
pub struct DmaBuilder {
    hull_length: Option<usize>,
    ema_length: Option<usize>,
    ema_gain_limit: Option<usize>,
    hull_ma_type: Option<String>,
    kernel: Kernel,
}

impl Default for DmaBuilder {
    fn default() -> Self {
        Self {
            hull_length: None,
            ema_length: None,
            ema_gain_limit: None,
            hull_ma_type: None,
            kernel: Kernel::Auto,
        }
    }
}

impl DmaBuilder {
    #[inline(always)]
    pub fn new() -> Self {
        Self::default()
    }

    #[inline(always)]
    pub fn hull_length(mut self, val: usize) -> Self {
        self.hull_length = Some(val);
        self
    }

    #[inline(always)]
    pub fn ema_length(mut self, val: usize) -> Self {
        self.ema_length = Some(val);
        self
    }

    #[inline(always)]
    pub fn ema_gain_limit(mut self, val: usize) -> Self {
        self.ema_gain_limit = Some(val);
        self
    }

    #[inline(always)]
    pub fn hull_ma_type(mut self, val: String) -> Self {
        self.hull_ma_type = Some(val);
        self
    }

    #[inline(always)]
    pub fn kernel(mut self, k: Kernel) -> Self {
        self.kernel = k;
        self
    }

    #[inline(always)]
    pub fn apply(self, c: &Candles) -> Result<DmaOutput, DmaError> {
        let p = DmaParams {
            hull_length: self.hull_length,
            ema_length: self.ema_length,
            ema_gain_limit: self.ema_gain_limit,
            hull_ma_type: self.hull_ma_type,
        };
        let i = DmaInput::from_candles(c, "close", p);
        dma_with_kernel(&i, self.kernel)
    }

    #[inline(always)]
    pub fn apply_slice(self, d: &[f64]) -> Result<DmaOutput, DmaError> {
        let p = DmaParams {
            hull_length: self.hull_length,
            ema_length: self.ema_length,
            ema_gain_limit: self.ema_gain_limit,
            hull_ma_type: self.hull_ma_type,
        };
        let i = DmaInput::from_slice(d, p);
        dma_with_kernel(&i, self.kernel)
    }

    #[inline(always)]
    pub fn into_stream(self) -> Result<DmaStream, DmaError> {
        let p = DmaParams {
            hull_length: self.hull_length,
            ema_length: self.ema_length,
            ema_gain_limit: self.ema_gain_limit,
            hull_ma_type: self.hull_ma_type,
        };
        DmaStream::try_new(p)
    }
}

#[derive(Debug, Error)]
pub enum DmaError {
    #[error("dma: Input data slice is empty.")]
    EmptyInputData,

    #[error("dma: All values are NaN.")]
    AllValuesNaN,

    #[error("dma: Invalid period: period = {period}, data length = {data_len}")]
    InvalidPeriod { period: usize, data_len: usize },

    #[error("dma: Not enough valid data: needed = {needed}, valid = {valid}")]
    NotEnoughValidData { needed: usize, valid: usize },

    #[error("dma: Invalid Hull MA type: {value}. Must be 'WMA' or 'EMA'.")]
    InvalidHullMAType { value: String },

    #[error("dma: Output slice length mismatch: expected = {expected}, got = {got}")]
    OutputLengthMismatch { expected: usize, got: usize },

    #[error("dma: Invalid range expansion: start = {start}, end = {end}, step = {step}")]
    InvalidRange {
        start: usize,
        end: usize,
        step: usize,
    },

    #[error("dma: Invalid kernel for batch path: {0:?}")]
    InvalidKernelForBatch(Kernel),
}

#[inline(always)]
pub fn dma(input: &DmaInput) -> Result<DmaOutput, DmaError> {
    dma_with_kernel(input, Kernel::Auto)
}

#[inline(always)]
pub fn dma_with_kernel(input: &DmaInput, kernel: Kernel) -> Result<DmaOutput, DmaError> {
    let (data, hull_len, ema_len, ema_gain_limit, hull_ma_type, first, chosen) =
        dma_prepare(input, kernel)?;

    let sqrt_len = (hull_len as f64).sqrt().round() as usize;
    let warmup_end = first + hull_len.max(ema_len) + sqrt_len - 1;

    let mut out = alloc_with_nan_prefix(data.len(), warmup_end);
    dma_compute_into(
        data,
        hull_len,
        ema_len,
        ema_gain_limit,
        &hull_ma_type,
        first,
        chosen,
        &mut out,
    );
    Ok(DmaOutput { values: out })
}

#[cfg(not(all(target_arch = "wasm32", feature = "wasm")))]
#[inline(always)]
pub fn dma_into(input: &DmaInput, out: &mut [f64]) -> Result<(), DmaError> {
    let (data, hull_len, ema_len, ema_gain_limit, hull_ma_type, first, chosen) =
        dma_prepare(input, Kernel::Auto)?;

    if out.len() != data.len() {
        return Err(DmaError::OutputLengthMismatch {
            expected: data.len(),
            got: out.len(),
        });
    }

    let sqrt_len = (hull_len as f64).sqrt().round() as usize;
    let warmup_end = first + hull_len.max(ema_len) + sqrt_len - 1;
    let end = warmup_end.min(out.len());
    let qnan = f64::from_bits(0x7ff8_0000_0000_0000);
    for v in &mut out[..end] {
        *v = qnan;
    }

    dma_compute_into(
        data,
        hull_len,
        ema_len,
        ema_gain_limit,
        &hull_ma_type,
        first,
        chosen,
        out,
    );
    Ok(())
}

#[inline(always)]
pub fn dma_into_slice(dst: &mut [f64], input: &DmaInput, kern: Kernel) -> Result<(), DmaError> {
    let (data, hull_len, ema_len, ema_gain_limit, hull_ma_type, first, chosen) =
        dma_prepare(input, kern)?;

    if dst.len() != data.len() {
        return Err(DmaError::OutputLengthMismatch {
            expected: data.len(),
            got: dst.len(),
        });
    }

    dma_compute_into(
        data,
        hull_len,
        ema_len,
        ema_gain_limit,
        &hull_ma_type,
        first,
        chosen,
        dst,
    );

    let sqrt_len = (hull_len as f64).sqrt().round() as usize;
    let warmup_end = first + hull_len.max(ema_len) + sqrt_len - 1;
    let end = warmup_end.min(dst.len());
    for v in &mut dst[..end] {
        *v = f64::NAN;
    }
    Ok(())
}

#[inline(always)]
fn dma_prepare<'a>(
    input: &'a DmaInput,
    kernel: Kernel,
) -> Result<(&'a [f64], usize, usize, usize, &'a str, usize, Kernel), DmaError> {
    let data: &[f64] = input.as_ref();
    let len = data.len();
    if len == 0 {
        return Err(DmaError::EmptyInputData);
    }

    let first = data
        .iter()
        .position(|x| !x.is_nan())
        .ok_or(DmaError::AllValuesNaN)?;
    let hull_length = input.get_hull_length();
    let ema_length = input.get_ema_length();
    let ema_gain_limit = input.get_ema_gain_limit();
    let hull_ma_type = input.hull_ma_type_str();

    if hull_length == 0 || hull_length > len {
        return Err(DmaError::InvalidPeriod {
            period: hull_length,
            data_len: len,
        });
    }
    if ema_length == 0 || ema_length > len {
        return Err(DmaError::InvalidPeriod {
            period: ema_length,
            data_len: len,
        });
    }

    let sqrt_len = (hull_length as f64).sqrt().round() as usize;
    let needed = hull_length.max(ema_length) + sqrt_len;
    if len - first < needed {
        return Err(DmaError::NotEnoughValidData {
            needed,
            valid: len - first,
        });
    }
    if hull_ma_type != "WMA" && hull_ma_type != "EMA" {
        return Err(DmaError::InvalidHullMAType {
            value: hull_ma_type.to_string(),
        });
    }
    let chosen = match kernel {
        Kernel::Auto => Kernel::Scalar,
        k => k,
    };
    Ok((
        data,
        hull_length,
        ema_length,
        ema_gain_limit,
        hull_ma_type,
        first,
        chosen,
    ))
}

#[inline(always)]
fn dma_compute_into(
    data: &[f64],
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    first: usize,
    kernel: Kernel,
    out: &mut [f64],
) {
    unsafe {
        #[cfg(all(target_arch = "wasm32", target_feature = "simd128"))]
        {
            if matches!(kernel, Kernel::Scalar | Kernel::ScalarBatch) {
                dma_simd128(
                    data,
                    hull_length,
                    ema_length,
                    ema_gain_limit,
                    hull_ma_type,
                    first,
                    out,
                );
                return;
            }
        }

        match kernel {
            Kernel::Scalar | Kernel::ScalarBatch => dma_scalar(
                data,
                hull_length,
                ema_length,
                ema_gain_limit,
                hull_ma_type,
                first,
                out,
            ),
            #[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
            Kernel::Avx2 | Kernel::Avx2Batch => dma_avx2(
                data,
                hull_length,
                ema_length,
                ema_gain_limit,
                hull_ma_type,
                first,
                out,
            ),
            #[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
            Kernel::Avx512 | Kernel::Avx512Batch => dma_avx512(
                data,
                hull_length,
                ema_length,
                ema_gain_limit,
                hull_ma_type,
                first,
                out,
            ),
            #[cfg(not(all(feature = "nightly-avx", target_arch = "x86_64")))]
            Kernel::Avx2 | Kernel::Avx2Batch | Kernel::Avx512 | Kernel::Avx512Batch => dma_scalar(
                data,
                hull_length,
                ema_length,
                ema_gain_limit,
                hull_ma_type,
                first,
                out,
            ),
            _ => unreachable!(),
        }
    }
}

#[inline]
pub fn dma_scalar(
    data: &[f64],
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    first: usize,
    out: &mut [f64],
) {
    let n = data.len();
    if n == 0 {
        return;
    }

    let alpha_e = 2.0 / (ema_length as f64 + 1.0);
    let one_minus_alpha_e = 1.0 - alpha_e;
    let i0_e = first + ema_length.saturating_sub(1);
    let mut e0_prev = 0.0;
    let mut e0_init_done = false;
    let mut ec_prev = 0.0;
    let mut ec_init_done = false;

    let half = hull_length / 2;
    let sqrt_len = (hull_length as f64).sqrt().round() as usize;

    let mut hull_val = f64::NAN;

    let wsum = |p: usize| -> f64 { (p * (p + 1)) as f64 / 2.0 };
    let i0_half = first + half.saturating_sub(1);
    let i0_full = first + hull_length.saturating_sub(1);

    let mut a_half = 0.0;
    let mut s_half = 0.0;
    let mut half_ready = false;

    let mut a_full = 0.0;
    let mut s_full = 0.0;
    let mut full_ready = false;

    let mut diff_ring: Vec<f64> = Vec::with_capacity(sqrt_len.max(1));
    let mut diff_pos: usize = 0;
    let mut diff_filled = 0usize;

    let mut a_diff = 0.0;
    let mut s_diff = 0.0;
    let mut diff_wma_init_done = false;

    let alpha_sqrt = if sqrt_len > 0 {
        2.0 / (sqrt_len as f64 + 1.0)
    } else {
        0.0
    };
    let mut diff_ema = 0.0;
    let mut diff_ema_init_done = false;
    let mut diff_sum_seed = 0.0;

    let mut e_half_prev = 0.0;
    let mut e_half_init_done = false;
    let mut e_full_prev = 0.0;
    let mut e_full_init_done = false;
    let alpha_half = if half > 0 {
        2.0 / (half as f64 + 1.0)
    } else {
        0.0
    };
    let alpha_full = if hull_length > 0 {
        2.0 / (hull_length as f64 + 1.0)
    } else {
        0.0
    };

    let is_wma = hull_ma_type == "WMA";

    for i in first..n {
        let x = data[i];

        if !e0_init_done {
            if i >= i0_e {
                let start = i + 1 - ema_length;
                let mut sum = 0.0;
                for k in start..=i {
                    sum += data[k];
                }
                e0_prev = sum / ema_length as f64;
                e0_init_done = true;
            }
        } else {
            e0_prev = x.mul_add(alpha_e, one_minus_alpha_e * e0_prev);
        }

        let mut diff_now = f64::NAN;

        if is_wma {
            if half > 0 {
                if !half_ready {
                    if i >= i0_half {
                        let start = i + 1 - half;
                        let mut sum = 0.0;
                        let mut wsum_local = 0.0;
                        for (j, idx) in (start..=i).enumerate() {
                            let w = (j + 1) as f64;
                            let v = data[idx];
                            sum += v;
                            wsum_local += w * v;
                        }
                        a_half = sum;
                        s_half = wsum_local;
                        half_ready = true;
                    }
                } else {
                    let a_prev = a_half;
                    a_half = a_prev + x - data[i - half];
                    s_half = s_half + (half as f64) * x - a_prev;
                }
            }

            if hull_length > 0 {
                if !full_ready {
                    if i >= i0_full {
                        let start = i + 1 - hull_length;
                        let mut sum = 0.0;
                        let mut wsum_local = 0.0;
                        for (j, idx) in (start..=i).enumerate() {
                            let w = (j + 1) as f64;
                            let v = data[idx];
                            sum += v;
                            wsum_local += w * v;
                        }
                        a_full = sum;
                        s_full = wsum_local;
                        full_ready = true;
                    }
                } else {
                    let a_prev = a_full;
                    a_full = a_prev + x - data[i - hull_length];
                    s_full = s_full + (hull_length as f64) * x - a_prev;
                }
            }

            if half_ready && full_ready {
                let w_half = s_half / wsum(half).max(1.0);
                let w_full = s_full / wsum(hull_length).max(1.0);
                diff_now = 2.0 * w_half - w_full;
            }
        } else {
            if half > 0 {
                if !e_half_init_done {
                    if i >= i0_half {
                        let start = i + 1 - half;
                        let mut sum = 0.0;
                        for k in start..=i {
                            sum += data[k];
                        }
                        e_half_prev = sum / half as f64;
                        e_half_init_done = true;
                    }
                } else {
                    e_half_prev = x.mul_add(alpha_half, (1.0 - alpha_half) * e_half_prev);
                }
            }

            if hull_length > 0 {
                if !e_full_init_done {
                    if i >= i0_full {
                        let start = i + 1 - hull_length;
                        let mut sum = 0.0;
                        for k in start..=i {
                            sum += data[k];
                        }
                        e_full_prev = sum / hull_length as f64;
                        e_full_init_done = true;
                    }
                } else {
                    e_full_prev = x.mul_add(alpha_full, (1.0 - alpha_full) * e_full_prev);
                }
            }

            if e_half_init_done && e_full_init_done {
                diff_now = 2.0 * e_half_prev - e_full_prev;
            }
        }

        if diff_now.is_finite() && sqrt_len > 0 {
            if diff_filled < sqrt_len {
                diff_ring.push(diff_now);
                diff_sum_seed += diff_now;
                diff_filled += 1;

                if diff_filled == sqrt_len {
                    if is_wma {
                        a_diff = 0.0;
                        s_diff = 0.0;
                        for (j, &v) in diff_ring.iter().enumerate() {
                            let w = (j + 1) as f64;
                            a_diff += v;
                            s_diff += w * v;
                        }
                        diff_wma_init_done = true;
                        hull_val = s_diff / wsum(sqrt_len).max(1.0);
                    } else {
                        diff_ema = diff_sum_seed / sqrt_len as f64;
                        diff_ema_init_done = true;
                        hull_val = diff_ema;
                    }
                }
            } else {
                let old = diff_ring[diff_pos];
                diff_ring[diff_pos] = diff_now;
                diff_pos = (diff_pos + 1) % sqrt_len;

                if is_wma {
                    let a_prev = a_diff;
                    a_diff = a_prev + diff_now - old;
                    s_diff = s_diff + (sqrt_len as f64) * diff_now - a_prev;
                    hull_val = s_diff / wsum(sqrt_len).max(1.0);
                } else {
                    diff_ema = diff_now.mul_add(alpha_sqrt, (1.0 - alpha_sqrt) * diff_ema);
                    hull_val = diff_ema;
                }
            }
        }

        let mut ec_now = f64::NAN;
        if e0_init_done {
            if !ec_init_done {
                ec_prev = e0_prev;
                ec_init_done = true;
                ec_now = ec_prev;
            } else {
                let dx = x - ec_prev;
                let t = alpha_e * dx;
                let base = e0_prev.mul_add(alpha_e, one_minus_alpha_e * ec_prev);
                let r = x - base;

                let g_sel = if t == 0.0 {
                    0.0
                } else {
                    let limit_i = ema_gain_limit as i64;
                    let target = (r / t) * 10.0;
                    let mut i0 = target.floor() as i64;
                    if i0 < 0 {
                        i0 = 0;
                    } else if i0 > limit_i {
                        i0 = limit_i;
                    }
                    let i1 = if i0 < limit_i { i0 + 1 } else { i0 };
                    let g0 = (i0 as f64) * 0.1;
                    let g1 = (i1 as f64) * 0.1;
                    let e0 = (r - t * g0).abs();
                    let e1 = (r - t * g1).abs();
                    if e0 <= e1 {
                        g0
                    } else {
                        g1
                    }
                };

                ec_now = (e0_prev + g_sel * dx).mul_add(alpha_e, one_minus_alpha_e * ec_prev);
                ec_prev = ec_now;
            }
        }

        if hull_val.is_finite() && ec_now.is_finite() {
            out[i] = 0.5 * (hull_val + ec_now);
        }
    }
}

#[cfg(all(target_arch = "wasm32", target_feature = "simd128"))]
#[inline]
unsafe fn dma_simd128(
    data: &[f64],
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    first_val: usize,
    out: &mut [f64],
) {
    use core::arch::wasm32::*;
    dma_scalar(
        data,
        hull_length,
        ema_length,
        ema_gain_limit,
        hull_ma_type,
        first_val,
        out,
    );
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn hsum256d(v: __m256d) -> f64 {
    let mut buf = [0.0f64; 4];
    _mm256_storeu_pd(buf.as_mut_ptr(), v);
    buf[0] + buf[1] + buf[2] + buf[3]
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn hsum512d(v: __m512d) -> f64 {
    let mut buf = [0.0f64; 8];
    _mm512_storeu_pd(buf.as_mut_ptr(), v);
    buf.iter().sum()
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn vabs256d(x: __m256d) -> __m256d {
    let sign = _mm256_set1_pd(-0.0);
    _mm256_andnot_pd(sign, x)
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn vabs512d(x: __m512d) -> __m512d {
    let sign = _mm512_set1_epi64(i64::MIN as i64);
    let xi = _mm512_castpd_si512(x);
    let cleared = _mm512_andnot_si512(sign, xi);
    _mm512_castsi512_pd(cleared)
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn sum_unweighted_avx2(ptr: *const f64, len: usize) -> f64 {
    let mut i = 0usize;
    let mut acc = _mm256_setzero_pd();
    while i + 4 <= len {
        let v = _mm256_loadu_pd(ptr.add(i));
        acc = _mm256_add_pd(acc, v);
        i += 4;
    }
    let mut s = hsum256d(acc);
    while i < len {
        s += *ptr.add(i);
        i += 1;
    }
    s
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn sum_unweighted_avx512(ptr: *const f64, len: usize) -> f64 {
    let mut i = 0usize;
    let mut acc = _mm512_setzero_pd();
    while i + 8 <= len {
        let v = _mm512_loadu_pd(ptr.add(i));
        acc = _mm512_add_pd(acc, v);
        i += 8;
    }
    let mut s = hsum512d(acc);
    while i < len {
        s += *ptr.add(i);
        i += 1;
    }
    s
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn seed_wma_window_avx2(ptr: *const f64, len: usize) -> (f64, f64) {
    let mut i = 0usize;
    let mut acc_v = _mm256_setzero_pd();
    let mut acc_wv = _mm256_setzero_pd();
    let inc = _mm256_set_pd(3.0, 2.0, 1.0, 0.0);
    let mut wbase = 1.0f64;
    while i + 4 <= len {
        let v = _mm256_loadu_pd(ptr.add(i));
        let w = _mm256_add_pd(_mm256_set1_pd(wbase), inc);
        acc_v = _mm256_add_pd(acc_v, v);
        acc_wv = _mm256_add_pd(acc_wv, _mm256_mul_pd(w, v));
        wbase += 4.0;
        i += 4;
    }
    let mut s = hsum256d(acc_v);
    let mut sw = hsum256d(acc_wv);
    while i < len {
        let val = *ptr.add(i);
        s += val;
        sw += (i as f64 + 1.0) * val;
        i += 1;
    }
    (s, sw)
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn seed_wma_window_avx512(ptr: *const f64, len: usize) -> (f64, f64) {
    let mut i = 0usize;
    let mut acc_v = _mm512_setzero_pd();
    let mut acc_wv = _mm512_setzero_pd();
    let inc = _mm512_set_pd(7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0);
    let mut wbase = 1.0f64;
    while i + 8 <= len {
        let v = _mm512_loadu_pd(ptr.add(i));
        let w = _mm512_add_pd(_mm512_set1_pd(wbase), inc);
        acc_v = _mm512_add_pd(acc_v, v);
        acc_wv = _mm512_add_pd(acc_wv, _mm512_mul_pd(w, v));
        wbase += 8.0;
        i += 8;
    }
    let mut s = hsum512d(acc_v);
    let mut sw = hsum512d(acc_wv);
    while i < len {
        let val = *ptr.add(i);
        s += val;
        sw += (i as f64 + 1.0) * val;
        i += 1;
    }
    (s, sw)
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn best_gain_search_avx2(
    x: f64,
    e0_prev: f64,
    ec_prev: f64,
    alpha_e: f64,
    ema_gain_limit: usize,
) -> f64 {
    let width = 4usize;
    let dx = _mm256_set1_pd(x - ec_prev);
    let x_v = _mm256_set1_pd(x);
    let e0_v = _mm256_set1_pd(e0_prev);
    let ec_prev_v = _mm256_set1_pd(ec_prev);
    let a_v = _mm256_set1_pd(alpha_e);
    let om_a_v = _mm256_set1_pd(1.0 - alpha_e);
    let inf_v = _mm256_set1_pd(f64::INFINITY);
    let limit_f = ema_gain_limit as f64;
    let limit_v = _mm256_set1_pd(limit_f);
    let scale = _mm256_set1_pd(0.1);

    let mut best_err = _mm256_set1_pd(f64::INFINITY);
    let mut best_g = _mm256_set1_pd(0.0);

    let mut idx = 0usize;
    while idx <= ema_gain_limit {
        let base = _mm256_set1_pd(idx as f64);
        let inc = _mm256_set_pd(3.0, 2.0, 1.0, 0.0);
        let idx_v = _mm256_add_pd(base, inc);

        let gt_mask = _mm256_cmp_pd(idx_v, limit_v, _CMP_GT_OQ);

        let g = _mm256_mul_pd(idx_v, scale);
        let e0_plus = _mm256_fmadd_pd(g, dx, e0_v);
        let pred = _mm256_fmadd_pd(a_v, e0_plus, _mm256_mul_pd(om_a_v, ec_prev_v));
        let err = vabs256d(_mm256_sub_pd(x_v, pred));

        let err_masked = _mm256_blendv_pd(err, inf_v, gt_mask);

        let lt = _mm256_cmp_pd(err_masked, best_err, _CMP_LT_OQ);
        best_err = _mm256_blendv_pd(best_err, err_masked, lt);
        best_g = _mm256_blendv_pd(best_g, g, lt);

        idx += width;
    }

    let mut e = [0.0f64; 4];
    let mut g = [0.0f64; 4];
    _mm256_storeu_pd(e.as_mut_ptr(), best_err);
    _mm256_storeu_pd(g.as_mut_ptr(), best_g);

    let mut best_e = f64::INFINITY;
    let mut best_gg = 0.0;
    for k in 0..4 {
        let ek = e[k];
        let gk = g[k];
        if ek < best_e || (ek == best_e && gk < best_gg) {
            best_e = ek;
            best_gg = gk;
        }
    }
    best_gg
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[inline(always)]
unsafe fn best_gain_search_avx512(
    x: f64,
    e0_prev: f64,
    ec_prev: f64,
    alpha_e: f64,
    ema_gain_limit: usize,
) -> f64 {
    let width = 8usize;
    let dx = _mm512_set1_pd(x - ec_prev);
    let x_v = _mm512_set1_pd(x);
    let e0_v = _mm512_set1_pd(e0_prev);
    let ec_prev_v = _mm512_set1_pd(ec_prev);
    let a_v = _mm512_set1_pd(alpha_e);
    let om_a_v = _mm512_set1_pd(1.0 - alpha_e);
    let inf_v = _mm512_set1_pd(f64::INFINITY);
    let limit_f = ema_gain_limit as f64;
    let limit_v = _mm512_set1_pd(limit_f);
    let scale = _mm512_set1_pd(0.1);

    let mut best_err = _mm512_set1_pd(f64::INFINITY);
    let mut best_g = _mm512_set1_pd(0.0);

    let mut idx = 0usize;
    while idx <= ema_gain_limit {
        let base = _mm512_set1_pd(idx as f64);
        let inc = _mm512_set_pd(7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0);
        let idx_v = _mm512_add_pd(base, inc);

        let k_invalid = _mm512_cmp_pd_mask(idx_v, limit_v, _CMP_GT_OQ);

        let g = _mm512_mul_pd(idx_v, scale);
        let e0_plus = _mm512_fmadd_pd(g, dx, e0_v);
        let pred = _mm512_fmadd_pd(a_v, e0_plus, _mm512_mul_pd(om_a_v, ec_prev_v));
        let err = vabs512d(_mm512_sub_pd(x_v, pred));

        let err_masked = _mm512_mask_mov_pd(err, k_invalid, inf_v);

        let k_lt = _mm512_cmp_pd_mask(err_masked, best_err, _CMP_LT_OQ);
        best_err = _mm512_mask_mov_pd(best_err, k_lt, err_masked);
        best_g = _mm512_mask_mov_pd(best_g, k_lt, g);

        idx += width;
    }

    let mut e = [0.0f64; 8];
    let mut g = [0.0f64; 8];
    _mm512_storeu_pd(e.as_mut_ptr(), best_err);
    _mm512_storeu_pd(g.as_mut_ptr(), best_g);

    let mut best_e = f64::INFINITY;
    let mut best_gg = 0.0;
    for k in 0..8 {
        let ek = e[k];
        let gk = g[k];
        if ek < best_e || (ek == best_e && gk < best_gg) {
            best_e = ek;
            best_gg = gk;
        }
    }
    best_gg
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[target_feature(enable = "avx2,fma")]
unsafe fn dma_avx2(
    data: &[f64],
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    first: usize,
    out: &mut [f64],
) {
    let n = data.len();
    if n == 0 {
        return;
    }

    let alpha_e = 2.0 / (ema_length as f64 + 1.0);
    let i0_e = first + ema_length.saturating_sub(1);
    let mut e0_prev = 0.0;
    let mut e0_init_done = false;
    let mut ec_prev = 0.0;
    let mut ec_init_done = false;

    let half = hull_length / 2;
    let sqrt_len = (hull_length as f64).sqrt().round() as usize;

    let mut hull_val = f64::NAN;

    let wsum = |p: usize| -> f64 { (p * (p + 1)) as f64 / 2.0 };
    let i0_half = first + half.saturating_sub(1);
    let i0_full = first + hull_length.saturating_sub(1);

    let mut a_half = 0.0;
    let mut s_half = 0.0;
    let mut half_ready = false;

    let mut a_full = 0.0;
    let mut s_full = 0.0;
    let mut full_ready = false;

    let mut diff_ring: Vec<f64> = Vec::with_capacity(sqrt_len.max(1));
    let mut diff_pos: usize = 0;
    let mut diff_filled = 0usize;

    let mut a_diff = 0.0;
    let mut s_diff = 0.0;
    let mut diff_wma_init_done = false;

    let alpha_sqrt = if sqrt_len > 0 {
        2.0 / (sqrt_len as f64 + 1.0)
    } else {
        0.0
    };
    let mut diff_ema = 0.0;
    let mut diff_ema_init_done = false;
    let mut diff_sum_seed = 0.0;

    let mut e_half_prev = 0.0;
    let mut e_half_init_done = false;
    let mut e_full_prev = 0.0;
    let mut e_full_init_done = false;
    let alpha_half = if half > 0 {
        2.0 / (half as f64 + 1.0)
    } else {
        0.0
    };
    let alpha_full = if hull_length > 0 {
        2.0 / (hull_length as f64 + 1.0)
    } else {
        0.0
    };

    let is_wma = hull_ma_type == "WMA";

    for i in first..n {
        let x = data[i];

        if !e0_init_done {
            if i >= i0_e {
                let start = i + 1 - ema_length;
                let sum = sum_unweighted_avx2(data.as_ptr().add(start), ema_length);
                e0_prev = sum / ema_length as f64;
                e0_init_done = true;
            }
        } else {
            e0_prev = x.mul_add(alpha_e, (1.0 - alpha_e) * e0_prev);
        }

        let mut diff_now = f64::NAN;

        if is_wma {
            if half > 0 {
                if !half_ready {
                    if i >= i0_half {
                        let start = i + 1 - half;
                        let (sum, wsum_local) =
                            seed_wma_window_avx2(data.as_ptr().add(start), half);
                        a_half = sum;
                        s_half = wsum_local;
                        half_ready = true;
                    }
                } else {
                    let a_prev = a_half;
                    a_half = a_prev + x - data[i - half];
                    s_half = s_half + (half as f64) * x - a_prev;
                }
            }

            if hull_length > 0 {
                if !full_ready {
                    if i >= i0_full {
                        let start = i + 1 - hull_length;
                        let (sum, wsum_local) =
                            seed_wma_window_avx2(data.as_ptr().add(start), hull_length);
                        a_full = sum;
                        s_full = wsum_local;
                        full_ready = true;
                    }
                } else {
                    let a_prev = a_full;
                    a_full = a_prev + x - data[i - hull_length];
                    s_full = s_full + (hull_length as f64) * x - a_prev;
                }
            }

            if half_ready && full_ready {
                let w_half = s_half / wsum(half).max(1.0);
                let w_full = s_full / wsum(hull_length).max(1.0);
                diff_now = 2.0 * w_half - w_full;
            }
        } else {
            if half > 0 {
                if !e_half_init_done {
                    if i >= i0_half {
                        let start = i + 1 - half;
                        let sum = sum_unweighted_avx2(data.as_ptr().add(start), half);
                        e_half_prev = sum / half as f64;
                        e_half_init_done = true;
                    }
                } else {
                    e_half_prev = x.mul_add(alpha_half, (1.0 - alpha_half) * e_half_prev);
                }
            }

            if hull_length > 0 {
                if !e_full_init_done {
                    if i >= i0_full {
                        let start = i + 1 - hull_length;
                        let sum = sum_unweighted_avx2(data.as_ptr().add(start), hull_length);
                        e_full_prev = sum / hull_length as f64;
                        e_full_init_done = true;
                    }
                } else {
                    e_full_prev = x.mul_add(alpha_full, (1.0 - alpha_full) * e_full_prev);
                }
            }

            if e_half_init_done && e_full_init_done {
                diff_now = 2.0 * e_half_prev - e_full_prev;
            }
        }

        if diff_now.is_finite() && sqrt_len > 0 {
            if diff_filled < sqrt_len {
                diff_ring.push(diff_now);
                diff_sum_seed += diff_now;
                diff_filled += 1;

                if diff_filled == sqrt_len {
                    if is_wma {
                        let (a0, s0) = seed_wma_window_avx2(diff_ring.as_ptr(), sqrt_len);
                        a_diff = a0;
                        s_diff = s0;
                        diff_wma_init_done = true;
                        let wsum_d = (sqrt_len * (sqrt_len + 1)) as f64 / 2.0;
                        hull_val = s_diff / wsum_d.max(1.0);
                    } else {
                        diff_ema = diff_sum_seed / sqrt_len as f64;
                        diff_ema_init_done = true;
                        hull_val = diff_ema;
                    }
                }
            } else {
                let old = diff_ring[diff_pos];
                diff_ring[diff_pos] = diff_now;
                diff_pos = (diff_pos + 1) % sqrt_len;

                if is_wma {
                    let a_prev = a_diff;
                    a_diff = a_prev + diff_now - old;
                    s_diff = s_diff + (sqrt_len as f64) * diff_now - a_prev;
                    let wsum_d = (sqrt_len * (sqrt_len + 1)) as f64 / 2.0;
                    hull_val = s_diff / wsum_d.max(1.0);
                } else {
                    diff_ema = diff_now.mul_add(alpha_sqrt, (1.0 - alpha_sqrt) * diff_ema);
                    hull_val = diff_ema;
                }
            }
        }

        let mut ec_now = f64::NAN;
        if e0_init_done {
            if !ec_init_done {
                ec_prev = e0_prev;
                ec_init_done = true;
                ec_now = ec_prev;
            } else {
                let dx = x - ec_prev;
                let t = alpha_e * dx;
                let base = e0_prev.mul_add(alpha_e, (1.0 - alpha_e) * ec_prev);
                let r = x - base;

                let g_sel = if t == 0.0 {
                    0.0
                } else {
                    let limit_i = ema_gain_limit as i64;
                    let target = (r / t) * 10.0;
                    let mut i0 = target.floor() as i64;
                    if i0 < 0 {
                        i0 = 0;
                    } else if i0 > limit_i {
                        i0 = limit_i;
                    }
                    let i1 = if i0 < limit_i { i0 + 1 } else { i0 };
                    let g0 = (i0 as f64) * 0.1;
                    let g1 = (i1 as f64) * 0.1;
                    let e0 = (r - t * g0).abs();
                    let e1 = (r - t * g1).abs();
                    if e0 <= e1 {
                        g0
                    } else {
                        g1
                    }
                };

                ec_now = (e0_prev + g_sel * dx).mul_add(alpha_e, (1.0 - alpha_e) * ec_prev);
                ec_prev = ec_now;
            }
        }

        if hull_val.is_finite() && ec_now.is_finite() {
            out[i] = 0.5 * (hull_val + ec_now);
        }
    }
}

#[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
#[target_feature(enable = "avx512f,fma")]
unsafe fn dma_avx512(
    data: &[f64],
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    first: usize,
    out: &mut [f64],
) {
    let n = data.len();
    if n == 0 {
        return;
    }

    let alpha_e = 2.0 / (ema_length as f64 + 1.0);
    let i0_e = first + ema_length.saturating_sub(1);
    let mut e0_prev = 0.0;
    let mut e0_init_done = false;
    let mut ec_prev = 0.0;
    let mut ec_init_done = false;

    let half = hull_length / 2;
    let sqrt_len = (hull_length as f64).sqrt().round() as usize;

    let mut hull_val = f64::NAN;

    let wsum = |p: usize| -> f64 { (p * (p + 1)) as f64 / 2.0 };
    let i0_half = first + half.saturating_sub(1);
    let i0_full = first + hull_length.saturating_sub(1);

    let mut a_half = 0.0;
    let mut s_half = 0.0;
    let mut half_ready = false;

    let mut a_full = 0.0;
    let mut s_full = 0.0;
    let mut full_ready = false;

    let mut diff_ring: Vec<f64> = Vec::with_capacity(sqrt_len.max(1));
    let mut diff_pos: usize = 0;
    let mut diff_filled = 0usize;

    let mut a_diff = 0.0;
    let mut s_diff = 0.0;
    let mut diff_wma_init_done = false;

    let alpha_sqrt = if sqrt_len > 0 {
        2.0 / (sqrt_len as f64 + 1.0)
    } else {
        0.0
    };
    let mut diff_ema = 0.0;
    let mut diff_ema_init_done = false;
    let mut diff_sum_seed = 0.0;

    let mut e_half_prev = 0.0;
    let mut e_half_init_done = false;
    let mut e_full_prev = 0.0;
    let mut e_full_init_done = false;
    let alpha_half = if half > 0 {
        2.0 / (half as f64 + 1.0)
    } else {
        0.0
    };
    let alpha_full = if hull_length > 0 {
        2.0 / (hull_length as f64 + 1.0)
    } else {
        0.0
    };

    let is_wma = hull_ma_type == "WMA";

    for i in first..n {
        let x = data[i];

        if !e0_init_done {
            if i >= i0_e {
                let start = i + 1 - ema_length;
                let sum = sum_unweighted_avx512(data.as_ptr().add(start), ema_length);
                e0_prev = sum / ema_length as f64;
                e0_init_done = true;
            }
        } else {
            e0_prev = x.mul_add(alpha_e, (1.0 - alpha_e) * e0_prev);
        }

        let mut diff_now = f64::NAN;

        if is_wma {
            if half > 0 {
                if !half_ready {
                    if i >= i0_half {
                        let start = i + 1 - half;
                        let (sum, wsum_local) =
                            seed_wma_window_avx512(data.as_ptr().add(start), half);
                        a_half = sum;
                        s_half = wsum_local;
                        half_ready = true;
                    }
                } else {
                    let a_prev = a_half;
                    a_half = a_prev + x - data[i - half];
                    s_half = s_half + (half as f64) * x - a_prev;
                }
            }

            if hull_length > 0 {
                if !full_ready {
                    if i >= i0_full {
                        let start = i + 1 - hull_length;
                        let (sum, wsum_local) =
                            seed_wma_window_avx512(data.as_ptr().add(start), hull_length);
                        a_full = sum;
                        s_full = wsum_local;
                        full_ready = true;
                    }
                } else {
                    let a_prev = a_full;
                    a_full = a_prev + x - data[i - hull_length];
                    s_full = s_full + (hull_length as f64) * x - a_prev;
                }
            }

            if half_ready && full_ready {
                let w_half = s_half / wsum(half).max(1.0);
                let w_full = s_full / wsum(hull_length).max(1.0);
                diff_now = 2.0 * w_half - w_full;
            }
        } else {
            if half > 0 {
                if !e_half_init_done {
                    if i >= i0_half {
                        let start = i + 1 - half;
                        let sum = sum_unweighted_avx512(data.as_ptr().add(start), half);
                        e_half_prev = sum / half as f64;
                        e_half_init_done = true;
                    }
                } else {
                    e_half_prev = x.mul_add(alpha_half, (1.0 - alpha_half) * e_half_prev);
                }
            }

            if hull_length > 0 {
                if !e_full_init_done {
                    if i >= i0_full {
                        let start = i + 1 - hull_length;
                        let sum = sum_unweighted_avx512(data.as_ptr().add(start), hull_length);
                        e_full_prev = sum / hull_length as f64;
                        e_full_init_done = true;
                    }
                } else {
                    e_full_prev = x.mul_add(alpha_full, (1.0 - alpha_full) * e_full_prev);
                }
            }

            if e_half_init_done && e_full_init_done {
                diff_now = 2.0 * e_half_prev - e_full_prev;
            }
        }

        if diff_now.is_finite() && sqrt_len > 0 {
            if diff_filled < sqrt_len {
                diff_ring.push(diff_now);
                diff_sum_seed += diff_now;
                diff_filled += 1;

                if diff_filled == sqrt_len {
                    if is_wma {
                        let (a0, s0) = seed_wma_window_avx512(diff_ring.as_ptr(), sqrt_len);
                        a_diff = a0;
                        s_diff = s0;
                        diff_wma_init_done = true;
                        let wsum_d = (sqrt_len * (sqrt_len + 1)) as f64 / 2.0;
                        hull_val = s_diff / wsum_d.max(1.0);
                    } else {
                        diff_ema = diff_sum_seed / sqrt_len as f64;
                        diff_ema_init_done = true;
                        hull_val = diff_ema;
                    }
                }
            } else {
                let old = diff_ring[diff_pos];
                diff_ring[diff_pos] = diff_now;
                diff_pos = (diff_pos + 1) % sqrt_len;

                if is_wma {
                    let a_prev = a_diff;
                    a_diff = a_prev + diff_now - old;
                    s_diff = s_diff + (sqrt_len as f64) * diff_now - a_prev;
                    let wsum_d = (sqrt_len * (sqrt_len + 1)) as f64 / 2.0;
                    hull_val = s_diff / wsum_d.max(1.0);
                } else {
                    diff_ema = diff_now.mul_add(alpha_sqrt, (1.0 - alpha_sqrt) * diff_ema);
                    hull_val = diff_ema;
                }
            }
        }

        let mut ec_now = f64::NAN;
        if e0_init_done {
            if !ec_init_done {
                ec_prev = e0_prev;
                ec_init_done = true;
                ec_now = ec_prev;
            } else {
                let dx = x - ec_prev;
                let t = alpha_e * dx;
                let base = e0_prev.mul_add(alpha_e, (1.0 - alpha_e) * ec_prev);
                let r = x - base;

                let g_sel = if t == 0.0 {
                    0.0
                } else {
                    let limit_i = ema_gain_limit as i64;
                    let target = (r / t) * 10.0;
                    let mut i0 = target.floor() as i64;
                    if i0 < 0 {
                        i0 = 0;
                    } else if i0 > limit_i {
                        i0 = limit_i;
                    }
                    let i1 = if i0 < limit_i { i0 + 1 } else { i0 };
                    let g0 = (i0 as f64) * 0.1;
                    let g1 = (i1 as f64) * 0.1;
                    let e0 = (r - t * g0).abs();
                    let e1 = (r - t * g1).abs();
                    if e0 <= e1 {
                        g0
                    } else {
                        g1
                    }
                };

                ec_now = (e0_prev + g_sel * dx).mul_add(alpha_e, (1.0 - alpha_e) * ec_prev);
                ec_prev = ec_now;
            }
        }

        if hull_val.is_finite() && ec_now.is_finite() {
            out[i] = 0.5 * (hull_val + ec_now);
        }
    }
}
#[derive(Debug, Clone)]
pub struct DmaStream {
    ema_length: usize,
    ema_gain_limit: usize,
    hull_length: usize,
    half: usize,
    sqrt_len: usize,
    is_wma: bool,

    cap: usize,
    ring: Vec<f64>,
    head: usize,
    filled: usize,

    i: usize,
    seen_first: bool,

    alpha_e: f64,
    sum_e0: f64,
    e0_prev: f64,
    e0_ready: bool,

    ec_prev: f64,
    ec_ready: bool,

    sum_half: f64,
    sum_full: f64,
    s_half: f64,
    s_full: f64,
    half_ready: bool,
    full_ready: bool,

    alpha_half: f64,
    alpha_full: f64,
    e_half_prev: f64,
    e_full_prev: f64,
    e_half_ready: bool,
    e_full_ready: bool,

    a_diff: f64,
    s_diff: f64,
    diff_wma_ready: bool,

    alpha_sqrt: f64,
    diff_ema: f64,
    diff_ema_ready: bool,

    diff_ring: Vec<f64>,
    diff_head: usize,
    diff_filled: usize,
}

impl DmaStream {
    pub fn try_new(params: DmaParams) -> Result<Self, DmaError> {
        let hull_length = params.hull_length.unwrap_or(7);
        let ema_length = params.ema_length.unwrap_or(20);
        let ema_gain_limit = params.ema_gain_limit.unwrap_or(50);
        let hull_ma_type = params.hull_ma_type.unwrap_or_else(|| "WMA".to_string());
        if hull_length == 0 || ema_length == 0 {
            return Err(DmaError::InvalidPeriod {
                period: hull_length.max(ema_length),
                data_len: 0,
            });
        }
        if hull_ma_type != "WMA" && hull_ma_type != "EMA" {
            return Err(DmaError::InvalidHullMAType {
                value: hull_ma_type,
            });
        }

        let half = hull_length / 2;
        let sqrt_len = (hull_length as f64).sqrt().round() as usize;
        let cap = hull_length.max(ema_length).max(1);

        Ok(Self {
            ema_length,
            ema_gain_limit,
            hull_length,
            half,
            sqrt_len,
            is_wma: hull_ma_type == "WMA",

            cap,
            ring: vec![f64::NAN; cap],
            head: 0,
            filled: 0,
            i: 0,
            seen_first: false,

            alpha_e: 2.0 / (ema_length as f64 + 1.0),
            sum_e0: 0.0,
            e0_prev: 0.0,
            e0_ready: false,

            ec_prev: 0.0,
            ec_ready: false,

            sum_half: 0.0,
            sum_full: 0.0,
            s_half: 0.0,
            s_full: 0.0,
            half_ready: half == 0,
            full_ready: hull_length == 0,

            alpha_half: if half > 0 {
                2.0 / (half as f64 + 1.0)
            } else {
                0.0
            },
            alpha_full: 2.0 / (hull_length as f64 + 1.0),
            e_half_prev: 0.0,
            e_full_prev: 0.0,
            e_half_ready: half == 0,
            e_full_ready: hull_length == 0,

            a_diff: 0.0,
            s_diff: 0.0,
            diff_wma_ready: sqrt_len == 0,
            alpha_sqrt: if sqrt_len > 0 {
                2.0 / (sqrt_len as f64 + 1.0)
            } else {
                0.0
            },
            diff_ema: 0.0,
            diff_ema_ready: sqrt_len == 0,
            diff_ring: vec![f64::NAN; sqrt_len.max(1)],
            diff_head: 0,
            diff_filled: 0,
        })
    }

    #[inline]
    pub fn update(&mut self, x: f64) -> Option<f64> {
        if !self.seen_first {
            self.i += 1;
            if x.is_nan() {
                return None;
            }
            self.seen_first = true;
        }

        let old_head = self.head;
        self.ring[old_head] = x;
        self.head = (old_head + 1) % self.cap;
        if self.filled < self.cap {
            self.filled += 1;
        }

        #[inline(always)]
        fn kback(ring: &[f64], head: usize, cap: usize, k: usize) -> f64 {
            let idx = (head + cap - k % cap) % cap;
            ring[idx]
        }

        if self.filled < self.ema_length {
            if x.is_finite() {
                self.sum_e0 += x;
            }
        } else {
            let out_e = kback(&self.ring, self.head, self.cap, self.ema_length);
            if x.is_finite() {
                self.sum_e0 += x;
            }
            if out_e.is_finite() {
                self.sum_e0 -= out_e;
            }
            if !self.e0_ready {
                self.e0_prev = self.sum_e0 / self.ema_length as f64;
                self.e0_ready = true;
            } else {
                self.e0_prev = self.alpha_e * x + (1.0 - self.alpha_e) * self.e0_prev;
            }
        }

        let mut diff_now = f64::NAN;

        if self.is_wma {
            if self.half > 0 {
                if self.filled < self.half {
                    if x.is_finite() {
                        self.sum_half += x;
                    }
                } else {
                    let out_h = kback(&self.ring, self.head, self.cap, self.half);
                    if x.is_finite() {
                        self.sum_half += x;
                    }
                    if out_h.is_finite() {
                        self.sum_half -= out_h;
                    }
                    if !self.half_ready {
                        self.s_half = 0.0;
                        for j in 0..self.half {
                            let v = kback(&self.ring, self.head, self.cap, self.half - j);
                            self.s_half += (j as f64 + 1.0) * v;
                        }
                        self.half_ready = true;
                    } else {
                        let a_prev = self.sum_half + out_h - x;
                        self.s_half = self.s_half + (self.half as f64) * x - a_prev;
                    }
                }
            } else {
                self.half_ready = true;
            }

            if self.filled < self.hull_length {
                if x.is_finite() {
                    self.sum_full += x;
                }
            } else {
                let out_f = kback(&self.ring, self.head, self.cap, self.hull_length);
                if x.is_finite() {
                    self.sum_full += x;
                }
                if out_f.is_finite() {
                    self.sum_full -= out_f;
                }
                if !self.full_ready {
                    self.s_full = 0.0;
                    for j in 0..self.hull_length {
                        let v = kback(&self.ring, self.head, self.cap, self.hull_length - j);
                        self.s_full += (j as f64 + 1.0) * v;
                    }
                    self.full_ready = true;
                } else {
                    let a_prev = self.sum_full + out_f - x;
                    self.s_full = self.s_full + (self.hull_length as f64) * x - a_prev;
                }
            }

            if self.half_ready && self.full_ready && self.sqrt_len > 0 {
                let wsum = |p: usize| (p * (p + 1)) as f64 / 2.0;
                let w_half = self.s_half / wsum(self.half).max(1.0);
                let w_full = self.s_full / wsum(self.hull_length).max(1.0);
                diff_now = 2.0 * w_half - w_full;
            }
        } else {
            if self.half > 0 {
                if self.filled < self.half {
                } else if !self.e_half_ready {
                    let mut s = 0.0;
                    for j in 0..self.half {
                        s += kback(&self.ring, self.head, self.cap, self.half - j);
                    }
                    self.e_half_prev = s / self.half as f64;
                    self.e_half_ready = true;
                } else {
                    self.e_half_prev =
                        self.alpha_half * x + (1.0 - self.alpha_half) * self.e_half_prev;
                }
            } else {
                self.e_half_ready = true;
            }

            if self.filled < self.hull_length {
            } else if !self.e_full_ready {
                let mut s = 0.0;
                for j in 0..self.hull_length {
                    s += kback(&self.ring, self.head, self.cap, self.hull_length - j);
                }
                self.e_full_prev = s / self.hull_length as f64;
                self.e_full_ready = true;
            } else {
                self.e_full_prev = self.alpha_full * x + (1.0 - self.alpha_full) * self.e_full_prev;
            }

            if self.e_half_ready && self.e_full_ready && self.sqrt_len > 0 {
                diff_now = 2.0 * self.e_half_prev - self.e_full_prev;
            }
        }

        let mut hull_val = f64::NAN;
        if self.sqrt_len == 0 {
            if diff_now.is_finite() {
                hull_val = diff_now;
            }
        } else if diff_now.is_finite() {
            let old = self.diff_ring[self.diff_head];
            self.diff_ring[self.diff_head] = diff_now;
            self.diff_head = (self.diff_head + 1) % self.sqrt_len;
            if self.diff_filled < self.sqrt_len {
                self.diff_filled += 1;
            }

            if self.is_wma {
                if !self.diff_wma_ready && self.diff_filled == self.sqrt_len {
                    self.a_diff = 0.0;
                    self.s_diff = 0.0;
                    for j in 0..self.sqrt_len {
                        let v = self.diff_ring[(self.diff_head + j) % self.sqrt_len];
                        self.a_diff += v;
                        self.s_diff += (j as f64 + 1.0) * v;
                    }
                    self.diff_wma_ready = true;
                    let wsum = (self.sqrt_len * (self.sqrt_len + 1)) as f64 / 2.0;
                    hull_val = self.s_diff / wsum.max(1.0);
                } else if self.diff_wma_ready {
                    let wsum = (self.sqrt_len * (self.sqrt_len + 1)) as f64 / 2.0;
                    let a_prev = self.a_diff;
                    self.s_diff = self.s_diff + (self.sqrt_len as f64) * diff_now - a_prev;
                    self.a_diff = a_prev + diff_now - old;
                    hull_val = self.s_diff / wsum.max(1.0);
                }
            } else {
                if !self.diff_ema_ready && self.diff_filled == self.sqrt_len {
                    let mut s = 0.0;
                    for j in 0..self.sqrt_len {
                        s += self.diff_ring[j];
                    }
                    self.diff_ema = s / self.sqrt_len as f64;
                    self.diff_ema_ready = true;
                    hull_val = self.diff_ema;
                } else if self.diff_ema_ready {
                    self.diff_ema =
                        self.alpha_sqrt * diff_now + (1.0 - self.alpha_sqrt) * self.diff_ema;
                    hull_val = self.diff_ema;
                }
            }
        }

        let mut ec_now = f64::NAN;
        if self.e0_ready {
            if !self.ec_ready {
                self.ec_prev = self.e0_prev;
                self.ec_ready = true;
                ec_now = self.ec_prev;
            } else {
                let one_minus_alpha_e = 1.0 - self.alpha_e;
                let dx = x - self.ec_prev;
                let t = self.alpha_e * dx;
                let base = self
                    .alpha_e
                    .mul_add(self.e0_prev, one_minus_alpha_e * self.ec_prev);
                let r = x - base;

                let g_sel = if t == 0.0 {
                    0.0
                } else {
                    let target = (r / t) * 10.0;
                    let limit_i = self.ema_gain_limit as i64;
                    let mut i0 = target.floor() as i64;
                    if i0 < 0 {
                        i0 = 0;
                    } else if i0 > limit_i {
                        i0 = limit_i;
                    }
                    let i1 = if i0 < limit_i { i0 + 1 } else { i0 };
                    let g0 = (i0 as f64) * 0.1;
                    let g1 = (i1 as f64) * 0.1;
                    let e0 = (r - t * g0).abs();
                    let e1 = (r - t * g1).abs();
                    if e0 <= e1 {
                        g0
                    } else {
                        g1
                    }
                };

                let ec = self.alpha_e.mul_add(
                    self.e0_prev + g_sel * dx,
                    (1.0 - self.alpha_e) * self.ec_prev,
                );
                self.ec_prev = ec;
                ec_now = ec;
            }
        }

        self.i += 1;

        if hull_val.is_finite() && ec_now.is_finite() {
            Some(0.5 * (hull_val + ec_now))
        } else {
            None
        }
    }
}

#[derive(Clone, Debug)]
pub struct DmaBatchRange {
    pub hull_length: (usize, usize, usize),
    pub ema_length: (usize, usize, usize),
    pub ema_gain_limit: (usize, usize, usize),
    pub hull_ma_type: String,
}

impl Default for DmaBatchRange {
    fn default() -> Self {
        Self {
            hull_length: (7, 7, 0),
            ema_length: (20, 269, 1),
            ema_gain_limit: (50, 50, 0),
            hull_ma_type: "WMA".to_string(),
        }
    }
}

#[derive(Clone, Debug, Default)]
pub struct DmaBatchBuilder {
    range: DmaBatchRange,
    kernel: Kernel,
}

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

    pub fn kernel(mut self, k: Kernel) -> Self {
        self.kernel = k;
        self
    }

    #[inline]
    pub fn hull_length_range(mut self, start: usize, end: usize, step: usize) -> Self {
        self.range.hull_length = (start, end, step);
        self
    }

    #[inline]
    pub fn hull_length_static(mut self, val: usize) -> Self {
        self.range.hull_length = (val, val, 0);
        self
    }

    #[inline]
    pub fn ema_length_range(mut self, start: usize, end: usize, step: usize) -> Self {
        self.range.ema_length = (start, end, step);
        self
    }

    #[inline]
    pub fn ema_length_static(mut self, val: usize) -> Self {
        self.range.ema_length = (val, val, 0);
        self
    }

    #[inline]
    pub fn ema_gain_limit_range(mut self, start: usize, end: usize, step: usize) -> Self {
        self.range.ema_gain_limit = (start, end, step);
        self
    }

    #[inline]
    pub fn ema_gain_limit_static(mut self, val: usize) -> Self {
        self.range.ema_gain_limit = (val, val, 0);
        self
    }

    #[inline]
    pub fn hull_ma_type(mut self, val: String) -> Self {
        self.range.hull_ma_type = val;
        self
    }

    pub fn apply_slice(self, data: &[f64]) -> Result<DmaBatchOutput, DmaError> {
        dma_batch_with_kernel(data, &self.range, self.kernel)
    }

    pub fn with_default_slice(data: &[f64], k: Kernel) -> Result<DmaBatchOutput, DmaError> {
        DmaBatchBuilder::new().kernel(k).apply_slice(data)
    }

    pub fn apply_candles(self, c: &Candles, src: &str) -> Result<DmaBatchOutput, DmaError> {
        let slice = source_type(c, src);
        self.apply_slice(slice)
    }

    pub fn with_default_candles(c: &Candles) -> Result<DmaBatchOutput, DmaError> {
        DmaBatchBuilder::new()
            .kernel(Kernel::Auto)
            .apply_candles(c, "close")
    }
}

#[derive(Clone, Debug)]
pub struct DmaBatchOutput {
    pub values: Vec<f64>,
    pub combos: Vec<DmaParams>,
    pub rows: usize,
    pub cols: usize,
}

impl DmaBatchOutput {
    pub fn row_for_params(&self, p: &DmaParams) -> Option<usize> {
        self.combos.iter().position(|c| {
            c.hull_length.unwrap_or(7) == p.hull_length.unwrap_or(7)
                && c.ema_length.unwrap_or(20) == p.ema_length.unwrap_or(20)
                && c.ema_gain_limit.unwrap_or(50) == p.ema_gain_limit.unwrap_or(50)
                && c.hull_ma_type.as_ref().unwrap_or(&"WMA".to_string())
                    == p.hull_ma_type.as_ref().unwrap_or(&"WMA".to_string())
        })
    }

    pub fn values_for(&self, p: &DmaParams) -> Option<&[f64]> {
        self.row_for_params(p).map(|row| {
            let start = row * self.cols;
            &self.values[start..start + self.cols]
        })
    }
}

#[inline(always)]
fn expand_grid_dma(r: &DmaBatchRange) -> Vec<DmaParams> {
    fn axis_usize((start, end, step): (usize, usize, usize)) -> Vec<usize> {
        if step == 0 || start == end {
            return vec![start];
        }
        if start < end {
            return (start..=end).step_by(step).collect();
        }

        let mut v: Vec<usize> = (end..=start).step_by(step).collect();
        v.reverse();
        v
    }

    let hull_lengths = axis_usize(r.hull_length);
    let ema_lengths = axis_usize(r.ema_length);
    let ema_gain_limits = axis_usize(r.ema_gain_limit);

    let mut combos = Vec::new();
    for &h in &hull_lengths {
        for &e in &ema_lengths {
            for &g in &ema_gain_limits {
                combos.push(DmaParams {
                    hull_length: Some(h),
                    ema_length: Some(e),
                    ema_gain_limit: Some(g),
                    hull_ma_type: Some(r.hull_ma_type.clone()),
                });
            }
        }
    }
    combos
}

#[inline(always)]
pub fn dma_batch_slice(
    data: &[f64],
    sweep: &DmaBatchRange,
    kern: Kernel,
) -> Result<DmaBatchOutput, DmaError> {
    dma_batch_inner(data, sweep, kern, false)
}

#[inline(always)]
pub fn dma_batch_par_slice(
    data: &[f64],
    sweep: &DmaBatchRange,
    kern: Kernel,
) -> Result<DmaBatchOutput, DmaError> {
    dma_batch_inner(data, sweep, kern, true)
}

#[inline(always)]
fn dma_batch_inner(
    data: &[f64],
    sweep: &DmaBatchRange,
    kern: Kernel,
    parallel: bool,
) -> Result<DmaBatchOutput, DmaError> {
    let combos = expand_grid_dma(sweep);
    let cols = data.len();
    let rows = combos.len();
    if cols == 0 {
        return Err(DmaError::EmptyInputData);
    }
    if rows == 0 {
        return Err(DmaError::InvalidRange {
            start: sweep.hull_length.0,
            end: sweep.hull_length.1,
            step: sweep.hull_length.2,
        });
    }

    let _cap = rows.checked_mul(cols).ok_or(DmaError::InvalidRange {
        start: rows,
        end: cols,
        step: 0,
    })?;

    let mut buf_mu = make_uninit_matrix(rows, cols);

    let first = data
        .iter()
        .position(|x| !x.is_nan())
        .ok_or(DmaError::AllValuesNaN)?;
    let warm: Vec<usize> = combos
        .iter()
        .map(|c| {
            let h = c.hull_length.unwrap();
            let e = c.ema_length.unwrap();
            let sqrt_len = (h as f64).sqrt().round() as usize;
            first + h.max(e) + sqrt_len - 1
        })
        .collect();
    init_matrix_prefixes(&mut buf_mu, cols, &warm);

    let mut guard = core::mem::ManuallyDrop::new(buf_mu);
    let out: &mut [f64] =
        unsafe { core::slice::from_raw_parts_mut(guard.as_mut_ptr() as *mut f64, guard.len()) };

    dma_batch_inner_into(data, sweep, kern, parallel, out)?;

    let values = unsafe {
        Vec::from_raw_parts(
            guard.as_mut_ptr() as *mut f64,
            guard.len(),
            guard.capacity(),
        )
    };

    Ok(DmaBatchOutput {
        values,
        combos,
        rows,
        cols,
    })
}

pub fn dma_batch_with_kernel(
    data: &[f64],
    sweep: &DmaBatchRange,
    k: Kernel,
) -> Result<DmaBatchOutput, DmaError> {
    let kernel = match k {
        Kernel::Auto => detect_best_batch_kernel(),
        other if other.is_batch() => other,
        other => return Err(DmaError::InvalidKernelForBatch(other)),
    };

    let simd = match kernel {
        Kernel::Avx512Batch => Kernel::Avx512,
        Kernel::Avx2Batch => Kernel::Avx2,
        Kernel::ScalarBatch => Kernel::Scalar,
        _ => unreachable!(),
    };
    dma_batch_par_slice(data, sweep, simd)
}

#[inline(always)]
fn dma_batch_inner_into(
    data: &[f64],
    sweep: &DmaBatchRange,
    k: Kernel,
    parallel: bool,
    out: &mut [f64],
) -> Result<Vec<DmaParams>, DmaError> {
    let combos = expand_grid_dma(sweep);
    if combos.is_empty() {
        return Err(DmaError::InvalidRange {
            start: sweep.hull_length.0,
            end: sweep.hull_length.1,
            step: sweep.hull_length.2,
        });
    }

    let first = data
        .iter()
        .position(|x| !x.is_nan())
        .ok_or(DmaError::AllValuesNaN)?;
    let cols = data.len();

    let actual = match k {
        Kernel::Auto => detect_best_batch_kernel(),
        other => other,
    };
    let simd = match actual {
        Kernel::Avx512Batch => Kernel::Avx512,
        Kernel::Avx2Batch => Kernel::Avx2,
        Kernel::ScalarBatch => Kernel::Scalar,
        _ => actual,
    };

    let do_row = |row: usize, dst_mu: &mut [MaybeUninit<f64>]| {
        let dst = unsafe {
            core::slice::from_raw_parts_mut(dst_mu.as_mut_ptr() as *mut f64, dst_mu.len())
        };
        let prm = &combos[row];
        let hull_len = prm.hull_length.unwrap_or(7);
        let ema_len = prm.ema_length.unwrap_or(20);

        let sqrt_len = (hull_len as f64).sqrt().round() as usize;
        let warmup_end = first + hull_len.max(ema_len) + sqrt_len - 1;
        let warmup_end = warmup_end.min(dst.len());

        for i in 0..warmup_end {
            dst[i] = f64::NAN;
        }

        dma_compute_into(
            data,
            hull_len,
            ema_len,
            prm.ema_gain_limit.unwrap_or(50),
            prm.hull_ma_type.as_ref().unwrap_or(&"WMA".to_string()),
            first,
            simd,
            dst,
        );
    };

    let dst_mu = unsafe {
        std::slice::from_raw_parts_mut(out.as_mut_ptr() as *mut MaybeUninit<f64>, out.len())
    };

    if parallel {
        #[cfg(not(target_arch = "wasm32"))]
        dst_mu
            .par_chunks_mut(cols)
            .enumerate()
            .for_each(|(r, row)| do_row(r, row));
        #[cfg(target_arch = "wasm32")]
        for (r, row) in dst_mu.chunks_mut(cols).enumerate() {
            do_row(r, row);
        }
    } else {
        for (r, row) in dst_mu.chunks_mut(cols).enumerate() {
            do_row(r, row);
        }
    }

    Ok(combos)
}

#[cfg(feature = "python")]
#[pyfunction(name = "dma")]
#[pyo3(signature = (data, hull_length=7, ema_length=20, ema_gain_limit=50, hull_ma_type="WMA", kernel=None))]
pub fn dma_py<'py>(
    py: Python<'py>,
    data: PyReadonlyArray1<'py, f64>,
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    kernel: Option<&str>,
) -> PyResult<Bound<'py, PyArray1<f64>>> {
    let slice_in = data.as_slice()?;
    let kern = validate_kernel(kernel, false)?;
    let params = DmaParams {
        hull_length: Some(hull_length),
        ema_length: Some(ema_length),
        ema_gain_limit: Some(ema_gain_limit),
        hull_ma_type: Some(hull_ma_type.to_string()),
    };
    let input = DmaInput::from_slice(slice_in, params);

    let result_vec: Vec<f64> = py
        .allow_threads(|| dma_with_kernel(&input, kern).map(|o| o.values))
        .map_err(|e| PyValueError::new_err(e.to_string()))?;

    Ok(result_vec.into_pyarray(py))
}

#[cfg(feature = "python")]
#[pyclass(name = "DmaStream")]
pub struct DmaStreamPy {
    stream: DmaStream,
}

#[cfg(feature = "python")]
#[pymethods]
impl DmaStreamPy {
    #[new]
    fn new(
        hull_length: usize,
        ema_length: usize,
        ema_gain_limit: usize,
        hull_ma_type: &str,
    ) -> PyResult<Self> {
        let params = DmaParams {
            hull_length: Some(hull_length),
            ema_length: Some(ema_length),
            ema_gain_limit: Some(ema_gain_limit),
            hull_ma_type: Some(hull_ma_type.to_string()),
        };
        let stream =
            DmaStream::try_new(params).map_err(|e| PyValueError::new_err(e.to_string()))?;
        Ok(DmaStreamPy { stream })
    }

    fn update(&mut self, value: f64) -> Option<f64> {
        self.stream.update(value)
    }
}

#[cfg(feature = "python")]
#[pyfunction(name = "dma_batch")]
#[pyo3(signature = (data, hull_length_range, ema_length_range, ema_gain_limit_range, hull_ma_type="WMA", kernel=None))]
pub fn dma_batch_py<'py>(
    py: Python<'py>,
    data: PyReadonlyArray1<'py, f64>,
    hull_length_range: (usize, usize, usize),
    ema_length_range: (usize, usize, usize),
    ema_gain_limit_range: (usize, usize, usize),
    hull_ma_type: &str,
    kernel: Option<&str>,
) -> PyResult<Bound<'py, PyDict>> {
    use numpy::{IntoPyArray, PyArray1, PyArrayMethods};
    let slice_in = data.as_slice()?;

    let sweep = DmaBatchRange {
        hull_length: hull_length_range,
        ema_length: ema_length_range,
        ema_gain_limit: ema_gain_limit_range,
        hull_ma_type: hull_ma_type.to_string(),
    };

    let combos = expand_grid_dma(&sweep);
    let rows = combos.len();
    let cols = slice_in.len();
    let total = rows
        .checked_mul(cols)
        .ok_or_else(|| PyValueError::new_err("rows*cols overflow"))?;

    let out_arr = unsafe { PyArray1::<f64>::new(py, [total], false) };
    let slice_out = unsafe { out_arr.as_slice_mut()? };

    let kern = validate_kernel(kernel, true)?;
    let combos = py
        .allow_threads(|| {
            let kernel = match kern {
                Kernel::Auto => detect_best_batch_kernel(),
                k => k,
            };
            let simd = match kernel {
                Kernel::Avx512Batch => Kernel::Avx512,
                Kernel::Avx2Batch => Kernel::Avx2,
                Kernel::ScalarBatch => Kernel::Scalar,
                other => return Err(DmaError::InvalidKernelForBatch(other)),
            };
            dma_batch_inner_into(slice_in, &sweep, simd, true, slice_out)
        })
        .map_err(|e| PyValueError::new_err(e.to_string()))?;

    let dict = PyDict::new(py);
    dict.set_item("values", out_arr.reshape((rows, cols))?)?;
    dict.set_item(
        "hull_lengths",
        combos
            .iter()
            .map(|p| p.hull_length.unwrap() as u64)
            .collect::<Vec<_>>()
            .into_pyarray(py),
    )?;
    dict.set_item(
        "ema_lengths",
        combos
            .iter()
            .map(|p| p.ema_length.unwrap() as u64)
            .collect::<Vec<_>>()
            .into_pyarray(py),
    )?;
    dict.set_item(
        "ema_gain_limits",
        combos
            .iter()
            .map(|p| p.ema_gain_limit.unwrap() as u64)
            .collect::<Vec<_>>()
            .into_pyarray(py),
    )?;
    dict.set_item("hull_ma_type", hull_ma_type)?;

    dict.set_item(
        "hull_ma_types",
        combos
            .iter()
            .map(|p| p.hull_ma_type.as_deref().unwrap_or("WMA"))
            .collect::<Vec<_>>(),
    )?;

    Ok(dict)
}

#[cfg(all(feature = "python", feature = "cuda"))]
#[pyfunction(name = "dma_cuda_batch_dev")]
#[pyo3(signature = (data_f32, hull_length_range, ema_length_range, ema_gain_limit_range, hull_ma_type="WMA", device_id=0))]
pub fn dma_cuda_batch_dev_py(
    py: Python<'_>,
    data_f32: numpy::PyReadonlyArray1<'_, f32>,
    hull_length_range: (usize, usize, usize),
    ema_length_range: (usize, usize, usize),
    ema_gain_limit_range: (usize, usize, usize),
    hull_ma_type: &str,
    device_id: usize,
) -> PyResult<DmaDeviceArrayF32Py> {
    if !cuda_available() {
        return Err(PyValueError::new_err("CUDA not available"));
    }

    let sweep = DmaBatchRange {
        hull_length: hull_length_range,
        ema_length: ema_length_range,
        ema_gain_limit: ema_gain_limit_range,
        hull_ma_type: hull_ma_type.to_string(),
    };

    let slice_in = data_f32.as_slice()?;
    let (inner, ctx, dev_id) = py.allow_threads(|| {
        let cuda = CudaDma::new(device_id).map_err(|e| PyValueError::new_err(e.to_string()))?;
        let ctx = cuda.context();
        let dev_id = cuda.device_id();
        let arr = cuda
            .dma_batch_dev(slice_in, &sweep)
            .map_err(|e| PyValueError::new_err(e.to_string()))?;
        Ok::<_, pyo3::PyErr>((arr, ctx, dev_id))
    })?;

    Ok(DmaDeviceArrayF32Py {
        inner: Some(inner),
        _ctx: ctx,
        device_id: dev_id,
    })
}

#[cfg(all(feature = "python", feature = "cuda"))]
#[pyfunction(name = "dma_cuda_many_series_one_param_dev")]
#[pyo3(signature = (data_tm_f32, hull_length, ema_length, ema_gain_limit, hull_ma_type="WMA", device_id=0))]
pub fn dma_cuda_many_series_one_param_dev_py(
    py: Python<'_>,
    data_tm_f32: PyReadonlyArray2<'_, f32>,
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
    device_id: usize,
) -> PyResult<DmaDeviceArrayF32Py> {
    if !cuda_available() {
        return Err(PyValueError::new_err("CUDA not available"));
    }

    let flat_in: &[f32] = data_tm_f32.as_slice()?;
    let rows = data_tm_f32.shape()[0];
    let cols = data_tm_f32.shape()[1];
    let params = DmaParams {
        hull_length: Some(hull_length),
        ema_length: Some(ema_length),
        ema_gain_limit: Some(ema_gain_limit),
        hull_ma_type: Some(hull_ma_type.to_string()),
    };

    let (inner, ctx, dev_id) = py.allow_threads(|| {
        let cuda = CudaDma::new(device_id).map_err(|e| PyValueError::new_err(e.to_string()))?;
        let ctx = cuda.context();
        let dev_id = cuda.device_id();
        let arr = cuda
            .dma_many_series_one_param_time_major_dev(flat_in, cols, rows, &params)
            .map_err(|e| PyValueError::new_err(e.to_string()))?;
        Ok::<_, pyo3::PyErr>((arr, ctx, dev_id))
    })?;

    Ok(DmaDeviceArrayF32Py {
        inner: Some(inner),
        _ctx: ctx,
        device_id: dev_id,
    })
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[wasm_bindgen]
pub fn dma_js(
    data: &[f64],
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
) -> Result<Vec<f64>, JsValue> {
    let params = DmaParams {
        hull_length: Some(hull_length),
        ema_length: Some(ema_length),
        ema_gain_limit: Some(ema_gain_limit),
        hull_ma_type: Some(hull_ma_type.to_string()),
    };
    let input = DmaInput::from_slice(data, params);

    let mut output = vec![0.0; data.len()];
    dma_into_slice(&mut output, &input, Kernel::Auto)
        .map_err(|e| JsValue::from_str(&e.to_string()))?;

    Ok(output)
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[wasm_bindgen]
pub fn dma_alloc(len: usize) -> *mut f64 {
    let mut vec = Vec::<f64>::with_capacity(len);
    let ptr = vec.as_mut_ptr();
    std::mem::forget(vec);
    ptr
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[wasm_bindgen]
pub fn dma_free(ptr: *mut f64, len: usize) {
    unsafe {
        let _ = Vec::from_raw_parts(ptr, len, len);
    }
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[wasm_bindgen]
pub fn dma_into(
    in_ptr: *const f64,
    out_ptr: *mut f64,
    len: usize,
    hull_length: usize,
    ema_length: usize,
    ema_gain_limit: usize,
    hull_ma_type: &str,
) -> Result<(), JsValue> {
    if in_ptr.is_null() || out_ptr.is_null() {
        return Err(JsValue::from_str("null pointer passed to dma_into"));
    }

    unsafe {
        let data = std::slice::from_raw_parts(in_ptr, len);

        let params = DmaParams {
            hull_length: Some(hull_length),
            ema_length: Some(ema_length),
            ema_gain_limit: Some(ema_gain_limit),
            hull_ma_type: Some(hull_ma_type.to_string()),
        };
        let input = DmaInput::from_slice(data, params);

        if in_ptr == out_ptr {
            let mut temp = vec![0.0; len];
            dma_into_slice(&mut temp, &input, Kernel::Auto)
                .map_err(|e| JsValue::from_str(&e.to_string()))?;
            let out = std::slice::from_raw_parts_mut(out_ptr, len);
            out.copy_from_slice(&temp);
        } else {
            let out = std::slice::from_raw_parts_mut(out_ptr, len);
            dma_into_slice(out, &input, Kernel::Auto)
                .map_err(|e| JsValue::from_str(&e.to_string()))?;
        }

        Ok(())
    }
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[derive(Serialize, Deserialize)]
pub struct DmaBatchConfig {
    pub hull_length_range: (usize, usize, usize),
    pub ema_length_range: (usize, usize, usize),
    pub ema_gain_limit_range: (usize, usize, usize),
    pub hull_ma_type: String,
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[derive(Serialize, Deserialize)]
pub struct DmaBatchJsOutput {
    pub values: Vec<f64>,
    pub combos: Vec<DmaParams>,
    pub rows: usize,
    pub cols: usize,
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[wasm_bindgen(js_name = dma_batch)]
pub fn dma_batch_unified_js(data: &[f64], config: JsValue) -> Result<JsValue, JsValue> {
    let cfg: DmaBatchConfig = serde_wasm_bindgen::from_value(config)
        .map_err(|e| JsValue::from_str(&format!("Invalid config: {e}")))?;

    let sweep = DmaBatchRange {
        hull_length: cfg.hull_length_range,
        ema_length: cfg.ema_length_range,
        ema_gain_limit: cfg.ema_gain_limit_range,
        hull_ma_type: cfg.hull_ma_type,
    };

    let combos = expand_grid_dma(&sweep);
    let rows = combos.len();
    let cols = data.len();
    if rows == 0 {
        return Err(JsValue::from_str("no parameter combinations"));
    }

    let mut buf_mu = make_uninit_matrix(rows, cols);
    let first = data
        .iter()
        .position(|x| !x.is_nan())
        .ok_or_else(|| JsValue::from_str("All NaN"))?;
    let warm: Vec<usize> = combos
        .iter()
        .map(|c| {
            let h = c.hull_length.unwrap();
            let e = c.ema_length.unwrap();
            let sqrt_len = (h as f64).sqrt().round() as usize;
            first + h.max(e) + sqrt_len - 1
        })
        .collect();
    init_matrix_prefixes(&mut buf_mu, cols, &warm);

    let mut guard = core::mem::ManuallyDrop::new(buf_mu);
    let out: &mut [f64] =
        unsafe { core::slice::from_raw_parts_mut(guard.as_mut_ptr() as *mut f64, guard.len()) };

    dma_batch_inner_into(data, &sweep, detect_best_kernel(), false, out)
        .map_err(|e| JsValue::from_str(&e.to_string()))?;

    let values = unsafe {
        Vec::from_raw_parts(
            guard.as_mut_ptr() as *mut f64,
            guard.len(),
            guard.capacity(),
        )
    };
    let js = DmaBatchJsOutput {
        values,
        combos,
        rows,
        cols,
    };
    serde_wasm_bindgen::to_value(&js)
        .map_err(|e| JsValue::from_str(&format!("Serialization error: {e}")))
}

#[cfg(all(target_arch = "wasm32", feature = "wasm"))]
#[wasm_bindgen]
pub fn dma_batch_into(
    in_ptr: *const f64,
    out_ptr: *mut f64,
    len: usize,
    hull_start: usize,
    hull_end: usize,
    hull_step: usize,
    ema_start: usize,
    ema_end: usize,
    ema_step: usize,
    gain_start: usize,
    gain_end: usize,
    gain_step: usize,
    hull_ma_type: &str,
) -> Result<usize, JsValue> {
    if in_ptr.is_null() || out_ptr.is_null() {
        return Err(JsValue::from_str("null pointer passed to dma_batch_into"));
    }
    unsafe {
        let data = std::slice::from_raw_parts(in_ptr, len);
        let sweep = DmaBatchRange {
            hull_length: (hull_start, hull_end, hull_step),
            ema_length: (ema_start, ema_end, ema_step),
            ema_gain_limit: (gain_start, gain_end, gain_step),
            hull_ma_type: hull_ma_type.to_string(),
        };
        let combos = expand_grid_dma(&sweep);
        let rows = combos.len();
        let cols = len;

        let out_mu = std::slice::from_raw_parts_mut(out_ptr as *mut MaybeUninit<f64>, rows * cols);
        let first = data
            .iter()
            .position(|x| !x.is_nan())
            .ok_or_else(|| JsValue::from_str("All NaN"))?;
        let warm: Vec<usize> = combos
            .iter()
            .map(|c| {
                let h = c.hull_length.unwrap();
                let e = c.ema_length.unwrap();
                let sqrt_len = (h as f64).sqrt().round() as usize;
                first + h.max(e) + sqrt_len - 1
            })
            .collect();
        init_matrix_prefixes(out_mu, cols, &warm);

        let out = std::slice::from_raw_parts_mut(out_ptr, rows * cols);
        dma_batch_inner_into(data, &sweep, detect_best_kernel(), false, out)
            .map_err(|e| JsValue::from_str(&e.to_string()))?;

        Ok(rows)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::skip_if_unsupported;
    use crate::utilities::data_loader::read_candles_from_csv;
    use std::error::Error;

    fn check_dma_accuracy(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let file_path = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let candles = read_candles_from_csv(file_path)?;

        let input = DmaInput::from_candles(&candles, "close", DmaParams::default());
        let result = dma_with_kernel(&input, kernel)?;

        let expected_last_five = [
            59404.62489256,
            59326.48766951,
            59195.35128538,
            59153.22811529,
            58933.88503421,
        ];

        let start = result.values.len().saturating_sub(5);
        for (i, &val) in result.values[start..].iter().enumerate() {
            let diff = (val - expected_last_five[i]).abs();
            assert!(
                diff < 0.001,
                "[{}] DMA {:?} mismatch at idx {}: got {}, expected {}, diff {}",
                test_name,
                kernel,
                i,
                val,
                expected_last_five[i],
                diff
            );
        }
        Ok(())
    }

    fn check_dma_partial_params(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let file_path = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let candles = read_candles_from_csv(file_path)?;

        let default_params = DmaParams {
            hull_length: None,
            ema_length: None,
            ema_gain_limit: None,
            hull_ma_type: None,
        };
        let input = DmaInput::from_candles(&candles, "close", default_params);
        let output = dma_with_kernel(&input, kernel)?;
        assert_eq!(output.values.len(), candles.close.len());

        Ok(())
    }

    fn check_dma_default_candles(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let file_path = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let candles = read_candles_from_csv(file_path)?;

        let input = DmaInput::with_default_candles(&candles);
        match input.data {
            DmaData::Candles { source, .. } => assert_eq!(source, "close"),
            _ => panic!("Expected DmaData::Candles"),
        }
        let output = dma_with_kernel(&input, kernel)?;
        assert_eq!(output.values.len(), candles.close.len());

        Ok(())
    }

    fn check_dma_zero_period(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let input_data = [10.0, 20.0, 30.0];
        let params = DmaParams {
            hull_length: Some(0),
            ema_length: None,
            ema_gain_limit: None,
            hull_ma_type: None,
        };
        let input = DmaInput::from_slice(&input_data, params);
        let res = dma_with_kernel(&input, kernel);
        assert!(
            res.is_err(),
            "[{}] DMA should fail with zero period",
            test_name
        );
        Ok(())
    }

    fn check_dma_period_exceeds_length(
        test_name: &str,
        kernel: Kernel,
    ) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let data_small = [10.0, 20.0, 30.0];
        let params = DmaParams {
            hull_length: Some(10),
            ema_length: None,
            ema_gain_limit: None,
            hull_ma_type: None,
        };
        let input = DmaInput::from_slice(&data_small, params);
        let res = dma_with_kernel(&input, kernel);
        assert!(
            res.is_err(),
            "[{}] DMA should fail with period exceeding length",
            test_name
        );
        Ok(())
    }

    fn check_dma_very_small_dataset(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let single_point = [42.0];
        let params = DmaParams::default();
        let input = DmaInput::from_slice(&single_point, params);
        let res = dma_with_kernel(&input, kernel);
        assert!(
            res.is_err(),
            "[{}] DMA should fail with insufficient data",
            test_name
        );
        Ok(())
    }

    fn check_dma_empty_input(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let empty: [f64; 0] = [];
        let params = DmaParams::default();
        let input = DmaInput::from_slice(&empty, params);
        let res = dma_with_kernel(&input, kernel);
        assert!(
            res.is_err(),
            "[{}] DMA should fail with empty input",
            test_name
        );
        Ok(())
    }

    fn check_dma_all_nan(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let nan_data = [f64::NAN, f64::NAN, f64::NAN];
        let params = DmaParams::default();
        let input = DmaInput::from_slice(&nan_data, params);
        let res = dma_with_kernel(&input, kernel);
        assert!(
            res.is_err(),
            "[{}] DMA should fail with all NaN values",
            test_name
        );
        Ok(())
    }

    fn check_dma_invalid_hull_type(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let data = [10.0; 50];
        let params = DmaParams {
            hull_length: Some(7),
            ema_length: Some(20),
            ema_gain_limit: Some(50),
            hull_ma_type: Some("INVALID".to_string()),
        };
        let input = DmaInput::from_slice(&data, params);
        let res = dma_with_kernel(&input, kernel);
        assert!(
            res.is_err(),
            "[{}] DMA should fail with invalid hull_ma_type",
            test_name
        );
        Ok(())
    }

    macro_rules! generate_all_dma_tests {
        ($($test_fn:ident),*) => {
            paste::paste! {
                $(
                    #[test] fn [<$test_fn _scalar>]() -> Result<(), Box<dyn Error>> { $test_fn(stringify!([<$test_fn _scalar>]), Kernel::Scalar) }
                    #[test] fn [<$test_fn _auto>  ]() -> Result<(), Box<dyn Error>> { $test_fn(stringify!([<$test_fn _auto>]),   Kernel::Auto) }
                )*
                #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                $(
                    #[test] fn [<$test_fn _avx2>  ]() -> Result<(), Box<dyn Error>> { $test_fn(stringify!([<$test_fn _avx2>]),   Kernel::Avx2) }
                    #[test] fn [<$test_fn _avx512>]() -> Result<(), Box<dyn Error>> { $test_fn(stringify!([<$test_fn _avx512>]), Kernel::Avx512) }
                )*
            }
        }
    }

    generate_all_dma_tests!(
        check_dma_accuracy,
        check_dma_partial_params,
        check_dma_default_candles,
        check_dma_zero_period,
        check_dma_period_exceeds_length,
        check_dma_very_small_dataset,
        check_dma_empty_input,
        check_dma_all_nan,
        check_dma_invalid_hull_type
    );

    macro_rules! generate_dma_batch_tests {
        ($($fn_name:ident),*) => {
            paste::paste! {
                $(
                    #[test]
                    fn [<$fn_name _scalar_batch>]() -> Result<(), Box<dyn Error>> {
                        $fn_name(stringify!([<$fn_name _scalar_batch>]), Kernel::ScalarBatch)
                    }
                )*
                #[cfg(all(feature = "nightly-avx", target_arch = "x86_64"))]
                $(
                    #[test]
                    fn [<$fn_name _avx2_batch>]() -> Result<(), Box<dyn Error>> {
                        $fn_name(stringify!([<$fn_name _avx2_batch>]), Kernel::Avx2Batch)
                    }
                    #[test]
                    fn [<$fn_name _avx512_batch>]() -> Result<(), Box<dyn Error>> {
                        $fn_name(stringify!([<$fn_name _avx512_batch>]), Kernel::Avx512Batch)
                    }
                )*
            }
        };
    }

    generate_dma_batch_tests!(check_dma_batch_basic);

    macro_rules! gen_batch_tests {
        ($fn_name:ident) => {
            paste::paste! {
                #[test] fn [<$fn_name _scalar>]()      { let _ = $fn_name(stringify!([<$fn_name _scalar>]), Kernel::ScalarBatch); }
                #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                #[test] fn [<$fn_name _avx2>]()        { let _ = $fn_name(stringify!([<$fn_name _avx2>]), Kernel::Avx2Batch); }
                #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                #[test] fn [<$fn_name _avx512>]()      { let _ = $fn_name(stringify!([<$fn_name _avx512>]), Kernel::Avx512Batch); }
                #[test] fn [<$fn_name _auto_detect>]() { let _ = $fn_name(stringify!([<$fn_name _auto_detect>]), Kernel::Auto); }
            }
        };
    }

    gen_batch_tests!(check_batch_sweep);

    fn check_dma_reinput(test: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test);
        let file = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let c = read_candles_from_csv(file)?;

        let first = DmaInput::from_candles(&c, "close", DmaParams::default());
        let out1 = dma_with_kernel(&first, kernel)?.values;

        let second = DmaInput::from_slice(&out1, DmaParams::default());
        let out2 = dma_with_kernel(&second, kernel)?.values;

        assert_eq!(out2.len(), out1.len());
        Ok(())
    }

    fn check_dma_nan_handling(test: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test);
        let file = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let c = read_candles_from_csv(file)?;

        let p = DmaParams::default();
        let input = DmaInput::from_candles(&c, "close", p.clone());
        let out = dma_with_kernel(&input, kernel)?.values;

        let first = c.close.iter().position(|x| !x.is_nan()).unwrap_or(0);
        let sqrt_len = (p.hull_length.unwrap_or(7) as f64).sqrt().round() as usize;
        let warm =
            first + p.hull_length.unwrap_or(7).max(p.ema_length.unwrap_or(20)) + sqrt_len - 1;
        for (i, &v) in out.iter().enumerate().skip(warm.min(out.len())) {
            assert!(!v.is_nan(), "[{test}] unexpected NaN at {i}");
        }
        Ok(())
    }

    fn check_batch_default_row(test: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test);
        let file = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let c = read_candles_from_csv(file)?;

        let out = DmaBatchBuilder::new()
            .kernel(kernel)
            .apply_candles(&c, "close")?;
        let def = DmaParams::default();
        let row = out.values_for(&def).expect("default row missing");
        assert_eq!(row.len(), c.close.len());
        Ok(())
    }

    fn check_batch_sweep(test: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test);
        let file = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let c = read_candles_from_csv(file)?;
        let out = DmaBatchBuilder::new()
            .kernel(kernel)
            .hull_length_range(7, 18, 1)
            .ema_length_range(10, 15, 1)
            .ema_gain_limit_range(10, 20, 5)
            .apply_candles(&c, "close")?;
        let expected = 12 * 6 * 3;
        assert_eq!(out.combos.len(), expected);
        assert_eq!(out.rows, expected);
        assert_eq!(out.cols, c.close.len());
        Ok(())
    }

    fn check_dma_streaming(test: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test);
        let file = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let c = read_candles_from_csv(file)?;
        let p = DmaParams::default();

        let batch =
            dma_with_kernel(&DmaInput::from_candles(&c, "close", p.clone()), kernel)?.values;

        let mut s = DmaStream::try_new(p)?;
        let mut stream = Vec::with_capacity(c.close.len());
        for &x in &c.close {
            stream.push(s.update(x).unwrap_or(f64::NAN));
        }

        assert_eq!(batch.len(), stream.len());
        for (i, (&b, &t)) in batch.iter().zip(&stream).enumerate() {
            if b.is_nan() && t.is_nan() {
                continue;
            }
            assert!(
                (b - t).abs() < 1e-9,
                "[{test}] idx {i} diff {}",
                (b - t).abs()
            );
        }
        Ok(())
    }

    macro_rules! gen_added_dma_tests {
        ($($f:ident),*) => {
            paste::paste! {
                $(
                    #[test] fn [<$f _scalar>]() -> Result<(), Box<dyn Error>> {
                        $f(stringify!([<$f _scalar>]), Kernel::Scalar)
                    }
                    #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                    #[test] fn [<$f _avx2>]() -> Result<(), Box<dyn Error>> {
                        $f(stringify!([<$f _avx2>]), Kernel::Avx2)
                    }
                    #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                    #[test] fn [<$f _avx512>]() -> Result<(), Box<dyn Error>> {
                        $f(stringify!([<$f _avx512>]), Kernel::Avx512)
                    }
                )*
            }
        }
    }

    gen_added_dma_tests!(check_dma_reinput, check_dma_nan_handling);

    macro_rules! gen_batch_sweep_tests {
        ($($f:ident),*) => {
            paste::paste! {
                $(
                    #[test] fn [<$f _scalar_batch>]() -> Result<(), Box<dyn Error>> {
                        $f(stringify!([<$f _scalar_batch>]), Kernel::ScalarBatch)
                    }
                    #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                    #[test] fn [<$f _avx2_batch>]() -> Result<(), Box<dyn Error>> {
                        $f(stringify!([<$f _avx2_batch>]), Kernel::Avx2Batch)
                    }
                    #[cfg(all(feature="nightly-avx", target_arch="x86_64"))]
                    #[test] fn [<$f _avx512_batch>]() -> Result<(), Box<dyn Error>> {
                        $f(stringify!([<$f _avx512_batch>]), Kernel::Avx512Batch)
                    }
                )*
            }
        }
    }

    gen_batch_sweep_tests!(check_batch_default_row);

    #[cfg(all(target_arch = "wasm32", target_feature = "simd128"))]
    #[test]
    fn test_dma_simd128_correctness() {
        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
        let p = DmaParams::default();
        let input = DmaInput::from_slice(&data, p);
        let scalar = dma_with_kernel(&input, Kernel::Scalar).unwrap();
        let simd = dma_with_kernel(&input, Kernel::Scalar).unwrap();
        assert_eq!(scalar.values.len(), simd.values.len());
        for (a, b) in scalar.values.iter().zip(simd.values.iter()) {
            assert!((a - b).abs() < 1e-10);
        }
    }

    #[cfg(debug_assertions)]
    #[test]
    fn test_dma_no_poison_values() -> Result<(), Box<dyn Error>> {
        let file_path = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let candles = read_candles_from_csv(file_path)?;

        let input = DmaInput::from_candles(&candles, "close", DmaParams::default());
        let output = dma(&input)?;

        for &v in &output.values {
            if v.is_nan() {
                continue;
            }
            let b = v.to_bits();

            assert_ne!(
                b, 0x11111111_11111111,
                "Found poison value 0x11111111_11111111"
            );
            assert_ne!(
                b, 0x22222222_22222222,
                "Found poison value 0x22222222_22222222"
            );
            assert_ne!(
                b, 0x33333333_33333333,
                "Found poison value 0x33333333_33333333"
            );
            assert_ne!(
                b, 0xDEADBEEF_DEADBEEF,
                "Found poison value 0xDEADBEEF_DEADBEEF"
            );
            assert_ne!(
                b, 0xFEEEFEEE_FEEEFEEE,
                "Found poison value 0xFEEEFEEE_FEEEFEEE"
            );
        }
        Ok(())
    }

    fn check_dma_batch_basic(test_name: &str, kernel: Kernel) -> Result<(), Box<dyn Error>> {
        skip_if_unsupported!(kernel, test_name);
        let data = vec![10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0];

        let sweep = DmaBatchRange {
            hull_length: (3, 5, 1),
            ema_length: (5, 5, 0),
            ema_gain_limit: (10, 10, 0),
            hull_ma_type: "WMA".to_string(),
        };
        let output = dma_batch_with_kernel(&data, &sweep, kernel)?;

        assert_eq!(
            output.rows, 3,
            "[{}] Expected 3 rows for hull_length range 3-5",
            test_name
        );
        assert_eq!(output.cols, data.len());
        assert_eq!(output.values.len(), output.rows * output.cols);
        assert_eq!(output.combos.len(), output.rows);

        Ok(())
    }

    #[test]
    fn test_dma_stream_incremental() -> Result<(), Box<dyn Error>> {
        let params = DmaParams {
            hull_length: Some(3),
            ema_length: Some(3),
            ema_gain_limit: Some(10),
            hull_ma_type: Some("WMA".to_string()),
        };

        let mut stream = DmaStream::try_new(params)?;
        let data = vec![10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0];

        let mut results = Vec::new();
        for &val in &data {
            if let Some(result) = stream.update(val) {
                results.push(result);
            }
        }

        assert!(
            !results.is_empty(),
            "Stream should produce results after warmup"
        );

        Ok(())
    }

    #[cfg(debug_assertions)]
    #[test]
    fn test_dma_batch_no_poison_values() -> Result<(), Box<dyn std::error::Error>> {
        use crate::utilities::data_loader::read_candles_from_csv;
        let file = "src/data/2018-09-01-2024-Bitfinex_Spot-4h.csv";
        let c = read_candles_from_csv(file)?;
        let out = DmaBatchBuilder::new()
            .hull_length_range(3, 8, 1)
            .ema_length_range(5, 10, 1)
            .ema_gain_limit_static(10)
            .apply_slice(&c.close)?;
        for &v in &out.values {
            if v.is_nan() {
                continue;
            }
            let b = v.to_bits();
            assert_ne!(b, 0x11111111_11111111);
            assert_ne!(b, 0x22222222_22222222);
            assert_ne!(b, 0x33333333_33333333);
        }
        Ok(())
    }

    #[test]
    fn test_dma_into_matches_api() -> Result<(), Box<dyn Error>> {
        let mut data = Vec::with_capacity(160);
        data.extend_from_slice(&[f64::NAN, f64::NAN, f64::NAN]);
        for i in 0..157 {
            let x = (i as f64 * 0.15).sin() * 5.0 + (i as f64) * 0.01;
            data.push(x);
        }

        let input = DmaInput::from_slice(&data, DmaParams::default());

        let baseline = dma(&input)?;

        let mut out = vec![0.0; data.len()];
        #[cfg(not(all(target_arch = "wasm32", feature = "wasm")))]
        {
            dma_into(&input, &mut out)?;
        }
        #[cfg(all(target_arch = "wasm32", feature = "wasm"))]
        {
            dma_into_slice(&mut out, &input, Kernel::Auto)?;
        }

        assert_eq!(baseline.values.len(), out.len());

        for (a, b) in baseline.values.iter().copied().zip(out.iter().copied()) {
            let both_nan = a.is_nan() && b.is_nan();
            assert!(both_nan || a == b, "mismatch: got {b:?}, expected {a:?}");
        }
        Ok(())
    }
}

#[cfg(all(feature = "python", feature = "cuda"))]
#[pyclass(module = "ta_indicators.cuda", name = "DmaDeviceArrayF32", unsendable)]
pub struct DmaDeviceArrayF32Py {
    pub(crate) inner: Option<DeviceArrayF32>,
    pub(crate) _ctx: Arc<Context>,
    pub(crate) device_id: u32,
}

#[cfg(all(feature = "python", feature = "cuda"))]
#[pymethods]
impl DmaDeviceArrayF32Py {
    #[getter]
    fn __cuda_array_interface__<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyDict>> {
        let inner = self
            .inner
            .as_ref()
            .ok_or_else(|| PyValueError::new_err("buffer already exported via __dlpack__"))?;
        let d = PyDict::new(py);
        d.set_item("shape", (inner.rows, inner.cols))?;
        d.set_item("typestr", "<f4")?;

        d.set_item(
            "strides",
            (
                inner.cols * std::mem::size_of::<f32>(),
                std::mem::size_of::<f32>(),
            ),
        )?;
        d.set_item("data", (inner.device_ptr() as usize, false))?;

        d.set_item("version", 3)?;
        Ok(d)
    }
    fn __dlpack_device__(&self) -> (i32, i32) {
        (2, self.device_id as i32)
    }

    #[pyo3(signature=(stream=None, max_version=None, dl_device=None, copy=None))]
    fn __dlpack__<'py>(
        &mut self,
        py: Python<'py>,
        stream: Option<pyo3::PyObject>,
        max_version: Option<pyo3::PyObject>,
        dl_device: Option<pyo3::PyObject>,
        copy: Option<pyo3::PyObject>,
    ) -> PyResult<pyo3::PyObject> {
        use crate::utilities::dlpack_cuda::export_f32_cuda_dlpack_2d;

        let (kdl, alloc_dev) = self.__dlpack_device__();
        if let Some(dev_obj) = dl_device.as_ref() {
            if let Ok((dev_ty, dev_id)) = dev_obj.extract::<(i32, i32)>(py) {
                if dev_ty != kdl || dev_id != alloc_dev {
                    let wants_copy = copy
                        .as_ref()
                        .and_then(|c| c.extract::<bool>(py).ok())
                        .unwrap_or(false);
                    if wants_copy {
                        return Err(PyValueError::new_err(
                            "device copy not implemented for __dlpack__",
                        ));
                    } else {
                        return Err(PyValueError::new_err("dl_device mismatch for __dlpack__"));
                    }
                }
            }
        }

        let _ = stream;

        let inner = self
            .inner
            .take()
            .ok_or_else(|| PyValueError::new_err("__dlpack__ may only be called once"))?;

        let rows = inner.rows;
        let cols = inner.cols;
        let buf = inner.buf;

        let max_version_bound = max_version.map(|obj| obj.into_bound(py));

        export_f32_cuda_dlpack_2d(py, buf, rows, cols, alloc_dev, max_version_bound)
    }
}