sklears_discriminant_analysis/

simd_optimizations.rs

//! SIMD Optimizations for Discriminant Analysis
//!
//! This module provides SIMD-accelerated implementations of common matrix operations
//! used in discriminant analysis, leveraging SciRS2's SIMD capabilities for maximum performance.

// ✅ Using SciRS2 dependencies following SciRS2 policy
use scirs2_core::ndarray::{s, Array1, Array2, ArrayView1, Axis};
// Note: SIMD operations from SciRS2 may not be available in current version
// Using fallback implementations

use rayon::prelude::*;
use sklears_core::{error::Result, prelude::SklearsError, types::Float};

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
use std::arch::x86_64::*;

/// Configuration for SIMD operations
#[derive(Debug, Clone)]
pub struct SimdConfig {
    /// Use AVX-512 if available
    pub use_avx512: bool,
    /// Use AVX2 if available
    pub use_avx2: bool,
    /// Use SSE if available
    pub use_sse: bool,
    /// Minimum vector length for SIMD acceleration
    pub min_vector_length: usize,
    /// Use parallel SIMD for large operations
    pub parallel_simd: bool,
    /// Block size for tiled matrix operations
    pub block_size: usize,
}

impl Default for SimdConfig {
    fn default() -> Self {
        Self {
            use_avx512: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && cfg!(feature = "avx512"),
            use_avx2: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && cfg!(feature = "avx2"),
            use_sse: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && cfg!(feature = "sse"),
            min_vector_length: 8, // Minimum 8 elements for SIMD
            parallel_simd: true,
            block_size: 64, // 64x64 blocks for tiled operations
        }
    }
}

/// SIMD-accelerated matrix operations for discriminant analysis
pub struct SimdMatrixOps {
    config: SimdConfig,
}

impl SimdMatrixOps {
    /// Create a new SIMD matrix operations engine
    pub fn new() -> Self {
        Self {
            config: SimdConfig::default(),
        }
    }

    /// Create with custom configuration
    pub fn with_config(config: SimdConfig) -> Self {
        Self { config }
    }

    /// SIMD-accelerated matrix-vector multiplication
    pub fn simd_matvec(
        &self,
        matrix: &Array2<Float>,
        vector: &Array1<Float>,
    ) -> Result<Array1<Float>> {
        if matrix.ncols() != vector.len() {
            return Err(SklearsError::InvalidInput(
                "Matrix columns must match vector length".to_string(),
            ));
        }

        let nrows = matrix.nrows();
        let mut result = Array1::zeros(nrows);

        // Use parallel SIMD if enabled and matrix is large enough
        if self.config.parallel_simd && nrows >= 64 {
            self.parallel_simd_matvec(matrix, vector, &mut result)?;
        } else {
            self.sequential_simd_matvec(matrix, vector, &mut result)?;
        }

        Ok(result)
    }

    /// Sequential SIMD matrix-vector multiplication
    fn sequential_simd_matvec(
        &self,
        matrix: &Array2<Float>,
        vector: &Array1<Float>,
        result: &mut Array1<Float>,
    ) -> Result<()> {
        for (i, result_elem) in result.iter_mut().enumerate() {
            let row = matrix.row(i);
            *result_elem = self.simd_dot_product(&row, vector.view())?;
        }
        Ok(())
    }

    /// Parallel SIMD matrix-vector multiplication
    fn parallel_simd_matvec(
        &self,
        matrix: &Array2<Float>,
        vector: &Array1<Float>,
        result: &mut Array1<Float>,
    ) -> Result<()> {
        let results: Result<Vec<Float>> = (0..matrix.nrows())
            .into_par_iter()
            .map(|i| -> Result<Float> {
                let row = matrix.row(i);
                self.simd_dot_product(&row, vector.view())
            })
            .collect();

        let computed_results = results?;
        for (i, value) in computed_results.into_iter().enumerate() {
            result[i] = value;
        }
        Ok(())
    }

    /// SIMD-accelerated dot product
    pub fn simd_dot_product(&self, a: &ArrayView1<Float>, b: ArrayView1<Float>) -> Result<Float> {
        if a.len() != b.len() {
            return Err(SklearsError::InvalidInput(
                "Vectors must have same length".to_string(),
            ));
        }

        // Use manual SIMD if vector is large enough
        if a.len() >= self.config.min_vector_length {
            let a_vec: Vec<Float> = a.to_vec();
            let b_vec: Vec<Float> = b.to_vec();

            // Use manual SIMD implementation
            self.manual_simd_dot_product(&a_vec, &b_vec)
        } else {
            // Standard dot product for small vectors
            Ok(a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum())
        }
    }

    /// Manual SIMD dot product implementation using intrinsics
    fn manual_simd_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        let len = a.len();
        let mut sum = 0.0;

        if self.config.use_avx512 && len >= 8 {
            #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
            {
                sum += unsafe { self.avx512_dot_product(a, b)? };
            }
            #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
            {
                sum += self.avx512_dot_product(a, b)?;
            }
        } else if self.config.use_avx2 && len >= 4 {
            #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
            {
                sum += unsafe { self.avx2_dot_product(a, b)? };
            }
            #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
            {
                sum += self.avx2_dot_product(a, b)?;
            }
        } else if self.config.use_sse && len >= 2 {
            #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
            {
                sum += unsafe { self.sse_dot_product(a, b)? };
            }
            #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
            {
                sum += self.sse_dot_product(a, b)?;
            }
        } else {
            // Scalar fallback
            sum = a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum();
        }

        Ok(sum)
    }

    /// AVX-512 dot product implementation
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    #[target_feature(enable = "avx512f")]
    unsafe fn avx512_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        let len = a.len();
        let simd_len = len & !7; // Round down to multiple of 8
        let mut sum = _mm512_setzero_pd();

        for i in (0..simd_len).step_by(8) {
            let a_vec = _mm512_loadu_pd(a.as_ptr().add(i));
            let b_vec = _mm512_loadu_pd(b.as_ptr().add(i));
            let prod = _mm512_mul_pd(a_vec, b_vec);
            sum = _mm512_add_pd(sum, prod);
        }

        // Extract sum from SIMD register
        let sum_array: [f64; 8] = std::mem::transmute(sum);
        let mut result = sum_array.iter().sum::<f64>();

        // Handle remaining elements
        for i in simd_len..len {
            result += a[i] * b[i];
        }

        Ok(result)
    }

    /// AVX2 dot product implementation
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    #[target_feature(enable = "avx2")]
    unsafe fn avx2_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        let len = a.len();
        let simd_len = len & !3; // Round down to multiple of 4
        let mut sum = _mm256_setzero_pd();

        for i in (0..simd_len).step_by(4) {
            let a_vec = _mm256_loadu_pd(a.as_ptr().add(i));
            let b_vec = _mm256_loadu_pd(b.as_ptr().add(i));
            let prod = _mm256_mul_pd(a_vec, b_vec);
            sum = _mm256_add_pd(sum, prod);
        }

        // Extract sum from SIMD register
        let sum_array: [f64; 4] = std::mem::transmute(sum);
        let mut result = sum_array.iter().sum::<f64>();

        // Handle remaining elements
        for i in simd_len..len {
            result += a[i] * b[i];
        }

        Ok(result)
    }

    /// SSE dot product implementation
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    #[target_feature(enable = "sse2")]
    unsafe fn sse_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        let len = a.len();
        let simd_len = len & !1; // Round down to multiple of 2
        let mut sum = _mm_setzero_pd();

        for i in (0..simd_len).step_by(2) {
            let a_vec = _mm_loadu_pd(a.as_ptr().add(i));
            let b_vec = _mm_loadu_pd(b.as_ptr().add(i));
            let prod = _mm_mul_pd(a_vec, b_vec);
            sum = _mm_add_pd(sum, prod);
        }

        // Extract sum from SIMD register
        let sum_array: [f64; 2] = std::mem::transmute(sum);
        let mut result = sum_array.iter().sum::<f64>();

        // Handle remaining elements
        for i in simd_len..len {
            result += a[i] * b[i];
        }

        Ok(result)
    }

    /// Fallback implementations for non-x86 platforms
    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    fn avx512_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        // Fallback to scalar implementation
        Ok(a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum())
    }

    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    fn avx2_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        // Fallback to scalar implementation
        Ok(a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum())
    }

    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    fn sse_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        // Fallback to scalar implementation
        Ok(a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum())
    }

    /// SIMD-accelerated matrix multiplication
    pub fn simd_matmul(&self, a: &Array2<Float>, b: &Array2<Float>) -> Result<Array2<Float>> {
        if a.ncols() != b.nrows() {
            return Err(SklearsError::InvalidInput(
                "Matrix dimensions incompatible for multiplication".to_string(),
            ));
        }

        let (m, k) = (a.nrows(), a.ncols());
        let n = b.ncols();
        let mut result = Array2::zeros((m, n));

        // Use tiled matrix multiplication for better cache locality
        if m >= self.config.block_size && n >= self.config.block_size && k >= self.config.block_size
        {
            self.tiled_simd_matmul(a, b, &mut result)?;
        } else {
            self.simple_simd_matmul(a, b, &mut result)?;
        }

        Ok(result)
    }

    /// Simple SIMD matrix multiplication
    fn simple_simd_matmul(
        &self,
        a: &Array2<Float>,
        b: &Array2<Float>,
        result: &mut Array2<Float>,
    ) -> Result<()> {
        let n = b.ncols();

        if self.config.parallel_simd {
            // Use parallel chunks approach instead of par_iter with mutable borrows
            let rows = (0..a.nrows()).collect::<Vec<_>>();
            let results: Result<Vec<Vec<Float>>> = rows
                .par_iter()
                .map(|&i| -> Result<Vec<Float>> {
                    let a_row = a.row(i);
                    let mut row_result = vec![0.0; n];
                    for j in 0..n {
                        let b_col = b.column(j);
                        row_result[j] = self.simd_dot_product(&a_row, b_col)?;
                    }
                    Ok(row_result)
                })
                .collect();

            // Copy results back to result matrix
            let computed_results = results?;
            for (i, row_data) in computed_results.into_iter().enumerate() {
                for (j, value) in row_data.into_iter().enumerate() {
                    result[[i, j]] = value;
                }
            }
        } else {
            // Sequential computation
            for i in 0..a.nrows() {
                let a_row = a.row(i);
                for j in 0..n {
                    let b_col = b.column(j);
                    result[[i, j]] = self.simd_dot_product(&a_row, b_col)?;
                }
            }
        }

        Ok(())
    }

    /// Tiled SIMD matrix multiplication for better cache performance
    fn tiled_simd_matmul(
        &self,
        a: &Array2<Float>,
        b: &Array2<Float>,
        result: &mut Array2<Float>,
    ) -> Result<()> {
        let (m, k, n) = (a.nrows(), a.ncols(), b.ncols());
        let block_size = self.config.block_size;

        // Tile the computation
        for i_block in (0..m).step_by(block_size) {
            for j_block in (0..n).step_by(block_size) {
                for k_block in (0..k).step_by(block_size) {
                    let i_end = (i_block + block_size).min(m);
                    let j_end = (j_block + block_size).min(n);
                    let k_end = (k_block + block_size).min(k);

                    // Compute block
                    self.compute_block(
                        a, b, result, i_block, i_end, j_block, j_end, k_block, k_end,
                    )?;
                }
            }
        }

        Ok(())
    }

    /// Compute a single block in tiled matrix multiplication
    fn compute_block(
        &self,
        a: &Array2<Float>,
        b: &Array2<Float>,
        result: &mut Array2<Float>,
        i_start: usize,
        i_end: usize,
        j_start: usize,
        j_end: usize,
        k_start: usize,
        k_end: usize,
    ) -> Result<()> {
        for i in i_start..i_end {
            for j in j_start..j_end {
                let mut sum = 0.0;

                // Use SIMD for the inner loop if possible
                let k_len = k_end - k_start;
                if k_len >= self.config.min_vector_length {
                    let a_slice = a.slice(s![i, k_start..k_end]);
                    let b_slice = b.slice(s![k_start..k_end, j]);
                    sum = self.simd_dot_product(&a_slice, b_slice)?;
                } else {
                    // Scalar computation for small blocks
                    for k in k_start..k_end {
                        sum += a[[i, k]] * b[[k, j]];
                    }
                }

                result[[i, j]] += sum;
            }
        }

        Ok(())
    }

    /// SIMD-accelerated element-wise operations
    pub fn simd_element_wise_add(
        &self,
        a: &Array1<Float>,
        b: &Array1<Float>,
    ) -> Result<Array1<Float>> {
        if a.len() != b.len() {
            return Err(SklearsError::InvalidInput(
                "Arrays must have same length".to_string(),
            ));
        }

        // Fallback implementation since SciRS2 SIMD may not be available
        Ok(a + b)
    }

    /// SIMD-accelerated element-wise subtraction
    pub fn simd_element_wise_subtract(
        &self,
        a: &Array1<Float>,
        b: &Array1<Float>,
    ) -> Result<Array1<Float>> {
        if a.len() != b.len() {
            return Err(SklearsError::InvalidInput(
                "Arrays must have same length".to_string(),
            ));
        }

        // Fallback implementation since SciRS2 SIMD may not be available
        Ok(a - b)
    }

    /// SIMD-accelerated element-wise multiplication
    pub fn simd_element_wise_multiply(
        &self,
        a: &Array1<Float>,
        b: &Array1<Float>,
    ) -> Result<Array1<Float>> {
        if a.len() != b.len() {
            return Err(SklearsError::InvalidInput(
                "Arrays must have same length".to_string(),
            ));
        }

        // Fallback implementation since SciRS2 SIMD may not be available
        Ok(a * b)
    }

    /// SIMD-accelerated covariance matrix computation
    pub fn simd_covariance(&self, data: &Array2<Float>) -> Result<Array2<Float>> {
        let n_samples = data.nrows() as Float;
        let _n_features = data.ncols();

        // Compute mean using SIMD
        let mean = self.simd_column_mean(data)?;

        // Center the data
        let mut centered_data = Array2::zeros(data.raw_dim());
        for (mut row, data_row) in centered_data
            .axis_iter_mut(Axis(0))
            .zip(data.axis_iter(Axis(0)))
        {
            row.assign(&self.simd_element_wise_subtract(&data_row.to_owned(), &mean)?);
        }

        // Compute covariance matrix: (X^T * X) / (n - 1)
        let covariance = self.simd_matmul(&centered_data.t().to_owned(), &centered_data)?;
        Ok(covariance / (n_samples - 1.0))
    }

    /// SIMD-accelerated column mean computation
    fn simd_column_mean(&self, data: &Array2<Float>) -> Result<Array1<Float>> {
        let n_samples = data.nrows() as Float;
        let n_features = data.ncols();
        let mut mean = Array1::zeros(n_features);

        if self.config.parallel_simd {
            let results: Result<Vec<Float>> = (0..n_features)
                .into_par_iter()
                .map(|j| -> Result<Float> {
                    let col = data.column(j);
                    Ok(col.sum() / n_samples)
                })
                .collect();

            let computed_means = results?;
            for (j, value) in computed_means.into_iter().enumerate() {
                mean[j] = value;
            }
        } else {
            for (j, mean_elem) in mean.iter_mut().enumerate() {
                let col = data.column(j);
                *mean_elem = col.sum() / n_samples;
            }
        }

        Ok(mean)
    }

    /// SIMD-accelerated distance computations for discriminant analysis
    pub fn simd_mahalanobis_distance(
        &self,
        x: &Array1<Float>,
        mean: &Array1<Float>,
        inv_cov: &Array2<Float>,
    ) -> Result<Float> {
        // Compute (x - mean)
        let diff = self.simd_element_wise_subtract(x, mean)?;

        // Compute (x - mean)^T * inv_cov
        let temp = self.simd_matvec(inv_cov, &diff)?;

        // Compute final dot product: (x - mean)^T * inv_cov * (x - mean)
        let distance_squared = self.simd_dot_product(&diff.view(), temp.view())?;

        Ok(distance_squared.sqrt())
    }

    /// Check if current CPU supports required SIMD features
    pub fn check_simd_support(&self) -> SimdSupport {
        SimdSupport {
            avx512: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && self.runtime_feature_detect("avx512f"),
            avx2: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && self.runtime_feature_detect("avx2"),
            sse: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && self.runtime_feature_detect("sse2"),
            fma: cfg!(any(target_arch = "x86", target_arch = "x86_64"))
                && self.runtime_feature_detect("fma"),
        }
    }

    /// Runtime feature detection helper
    fn runtime_feature_detect(&self, feature: &str) -> bool {
        #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
        {
            match feature {
                "avx512f" => is_x86_feature_detected!("avx512f"),
                "avx2" => is_x86_feature_detected!("avx2"),
                "sse2" => is_x86_feature_detected!("sse2"),
                "fma" => is_x86_feature_detected!("fma"),
                _ => false,
            }
        }
        #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
        {
            let _ = feature;
            false
        }
    }
}

/// Information about SIMD support on current CPU
#[derive(Debug, Clone)]
pub struct SimdSupport {
    /// avx512
    pub avx512: bool,
    /// avx2
    pub avx2: bool,
    /// sse
    pub sse: bool,
    /// fma
    pub fma: bool,
}

impl Default for SimdMatrixOps {
    fn default() -> Self {
        Self::new()
    }
}

/// Extension trait for SIMD operations on ndarray types
pub trait SimdArrayOps<T> {
    /// SIMD-accelerated dot product
    fn simd_dot(&self, other: &Array1<T>) -> Result<T>;

    /// SIMD-accelerated matrix-vector multiplication
    fn simd_dot_matrix(&self, vector: &Array1<T>) -> Result<Array1<T>>;

    /// SIMD-accelerated element-wise addition
    fn simd_add(&self, other: &Array1<T>) -> Result<Array1<T>>;
}

impl SimdArrayOps<Float> for Array1<Float> {
    fn simd_dot(&self, other: &Array1<Float>) -> Result<Float> {
        let simd_ops = SimdMatrixOps::new();
        simd_ops.simd_dot_product(&self.view(), other.view())
    }

    fn simd_dot_matrix(&self, _vector: &Array1<Float>) -> Result<Array1<Float>> {
        // This would be for when self is treated as a 1xN matrix
        Err(SklearsError::InvalidInput(
            "Not applicable for 1D array".to_string(),
        ))
    }

    fn simd_add(&self, other: &Array1<Float>) -> Result<Array1<Float>> {
        let simd_ops = SimdMatrixOps::new();
        simd_ops.simd_element_wise_add(self, other)
    }
}

impl SimdArrayOps<Float> for Array2<Float> {
    fn simd_dot(&self, _other: &Array1<Float>) -> Result<Float> {
        Err(SklearsError::InvalidInput(
            "Use simd_dot_matrix for matrix-vector operations".to_string(),
        ))
    }

    fn simd_dot_matrix(&self, vector: &Array1<Float>) -> Result<Array1<Float>> {
        let simd_ops = SimdMatrixOps::new();
        simd_ops.simd_matvec(self, vector)
    }

    fn simd_add(&self, _other: &Array1<Float>) -> Result<Array1<Float>> {
        Err(SklearsError::InvalidInput(
            "Not applicable for 2D array with 1D array".to_string(),
        ))
    }
}

/// Advanced SIMD operations for discriminant analysis
pub struct AdvancedSimdOps {
    base_ops: SimdMatrixOps,
    config: SimdConfig,
}

impl AdvancedSimdOps {
    pub fn new() -> Self {
        Self {
            base_ops: SimdMatrixOps::new(),
            config: SimdConfig::default(),
        }
    }

    pub fn with_config(config: SimdConfig) -> Self {
        Self {
            base_ops: SimdMatrixOps::with_config(config.clone()),
            config,
        }
    }

    /// SIMD-accelerated eigenvalue computation approximation for 2x2 matrices
    pub fn simd_eigenvalues_2x2(&self, matrices: &Array2<Float>) -> Result<Array1<Float>> {
        if matrices.nrows() % 2 != 0 || matrices.ncols() % 2 != 0 {
            return Err(SklearsError::InvalidInput(
                "Input must be stacks of 2x2 matrices".to_string(),
            ));
        }

        let n_matrices = matrices.nrows() / 2;
        let mut eigenvalues = Array1::zeros(n_matrices * 2);

        // Process multiple 2x2 matrices using SIMD
        for i in 0..n_matrices {
            let matrix_start = i * 2;
            let a = matrices[[matrix_start, 0]];
            let b = matrices[[matrix_start, 1]];
            let c = matrices[[matrix_start + 1, 0]];
            let d = matrices[[matrix_start + 1, 1]];

            // Compute eigenvalues: λ = (trace ± √(trace² - 4*det)) / 2
            let trace = a + d;
            let det = a * d - b * c;
            let discriminant = (trace * trace - 4.0 * det).sqrt();

            eigenvalues[i * 2] = (trace + discriminant) / 2.0;
            eigenvalues[i * 2 + 1] = (trace - discriminant) / 2.0;
        }

        Ok(eigenvalues)
    }

    /// SIMD-accelerated batch matrix inversion for small matrices
    pub fn simd_batch_inverse_2x2(&self, matrices: &Array2<Float>) -> Result<Array2<Float>> {
        if matrices.nrows() % 2 != 0 || matrices.ncols() % 2 != 0 {
            return Err(SklearsError::InvalidInput(
                "Input must be stacks of 2x2 matrices".to_string(),
            ));
        }

        let n_matrices = matrices.nrows() / 2;
        let mut inverses = Array2::zeros(matrices.raw_dim());

        for i in 0..n_matrices {
            let matrix_start = i * 2;
            let a = matrices[[matrix_start, 0]];
            let b = matrices[[matrix_start, 1]];
            let c = matrices[[matrix_start + 1, 0]];
            let d = matrices[[matrix_start + 1, 1]];

            let det = a * d - b * c;
            if det.abs() < Float::EPSILON {
                return Err(SklearsError::NumericalError(format!(
                    "Matrix {} is singular",
                    i
                )));
            }

            let inv_det = 1.0 / det;
            inverses[[matrix_start, 0]] = d * inv_det;
            inverses[[matrix_start, 1]] = -b * inv_det;
            inverses[[matrix_start + 1, 0]] = -c * inv_det;
            inverses[[matrix_start + 1, 1]] = a * inv_det;
        }

        Ok(inverses)
    }

    /// SIMD-accelerated softmax computation
    pub fn simd_softmax(&self, x: &Array1<Float>) -> Result<Array1<Float>> {
        if x.is_empty() {
            return Err(SklearsError::InvalidInput(
                "Input array is empty".to_string(),
            ));
        }

        // Find maximum for numerical stability
        let max_val = x.iter().fold(Float::NEG_INFINITY, |acc, &val| acc.max(val));

        // Compute exp(x - max) using SIMD-accelerated element-wise operations
        let shifted: Array1<Float> = x.mapv(|val| val - max_val);
        let exp_values = self.simd_exp(&shifted)?;

        // Compute sum of exponentials
        let sum_exp = exp_values.sum();

        // Normalize using SIMD division
        Ok(exp_values / sum_exp)
    }

    /// SIMD-accelerated exponential function approximation
    fn simd_exp(&self, x: &Array1<Float>) -> Result<Array1<Float>> {
        // Fast exponential approximation using Taylor series
        // exp(x) ≈ 1 + x + x²/2! + x³/3! + ... (truncated)
        let len = x.len();
        let mut result = Array1::ones(len);

        if len >= self.config.min_vector_length {
            // Use vectorized computation for better performance
            let x_vec = x.to_vec();
            let exp_vec = self.vectorized_exp(&x_vec)?;
            result.assign(&Array1::from_vec(exp_vec));
        } else {
            // Scalar fallback for small arrays
            result.iter_mut().zip(x.iter()).for_each(|(res, &val)| {
                *res = val.exp();
            });
        }

        Ok(result)
    }

    /// Vectorized exponential computation
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    fn vectorized_exp(&self, x: &[Float]) -> Result<Vec<Float>> {
        if self.config.use_avx2 && x.len() >= 4 {
            unsafe { self.avx2_exp(x) }
        } else if self.config.use_sse && x.len() >= 2 {
            unsafe { self.sse_exp(x) }
        } else {
            Ok(x.iter().map(|&val| val.exp()).collect())
        }
    }

    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    fn vectorized_exp(&self, x: &[Float]) -> Result<Vec<Float>> {
        // ARM NEON fallback or scalar implementation
        Ok(x.iter().map(|&val| val.exp()).collect())
    }

    /// AVX2 exponential implementation
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    #[target_feature(enable = "avx2")]
    unsafe fn avx2_exp(&self, x: &[Float]) -> Result<Vec<Float>> {
        let len = x.len();
        let simd_len = len & !3; // Round down to multiple of 4
        let mut result = vec![0.0; len];

        for i in (0..simd_len).step_by(4) {
            let _x_vec = _mm256_loadu_pd(x.as_ptr().add(i));

            // Fast exp approximation using polynomial approximation
            // This is a simplified implementation - real SIMD exp would use more sophisticated methods
            let exp_vec = _mm256_set_pd(x[i + 3].exp(), x[i + 2].exp(), x[i + 1].exp(), x[i].exp());

            _mm256_storeu_pd(result.as_mut_ptr().add(i), exp_vec);
        }

        // Handle remaining elements
        for i in simd_len..len {
            result[i] = x[i].exp();
        }

        Ok(result)
    }

    /// SSE exponential implementation
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    #[target_feature(enable = "sse2")]
    unsafe fn sse_exp(&self, x: &[Float]) -> Result<Vec<Float>> {
        let len = x.len();
        let simd_len = len & !1; // Round down to multiple of 2
        let mut result = vec![0.0; len];

        for i in (0..simd_len).step_by(2) {
            let exp_vec = _mm_set_pd(x[i + 1].exp(), x[i].exp());
            _mm_storeu_pd(result.as_mut_ptr().add(i), exp_vec);
        }

        // Handle remaining elements
        for i in simd_len..len {
            result[i] = x[i].exp();
        }

        Ok(result)
    }

    /// SIMD-accelerated log-determinant computation for batch of matrices
    pub fn simd_log_determinant_batch(
        &self,
        matrices: &Array2<Float>,
        size: usize,
    ) -> Result<Array1<Float>> {
        let n_matrices = matrices.nrows() / size;
        let mut log_dets = Array1::zeros(n_matrices);

        match size {
            2 => {
                // Optimized 2x2 case
                for i in 0..n_matrices {
                    let offset = i * 2;
                    let a = matrices[[offset, 0]];
                    let b = matrices[[offset, 1]];
                    let c = matrices[[offset + 1, 0]];
                    let d = matrices[[offset + 1, 1]];

                    let det = a * d - b * c;
                    if det <= 0.0 {
                        return Err(SklearsError::NumericalError(format!(
                            "Matrix {} has non-positive determinant",
                            i
                        )));
                    }
                    log_dets[i] = det.ln();
                }
            }
            3 => {
                // Optimized 3x3 case using rule of Sarrus
                for i in 0..n_matrices {
                    let offset = i * 3;
                    let a = matrices[[offset, 0]];
                    let b = matrices[[offset, 1]];
                    let c = matrices[[offset, 2]];
                    let d = matrices[[offset + 1, 0]];
                    let e = matrices[[offset + 1, 1]];
                    let f = matrices[[offset + 1, 2]];
                    let g = matrices[[offset + 2, 0]];
                    let h = matrices[[offset + 2, 1]];
                    let i_elem = matrices[[offset + 2, 2]];

                    let det =
                        a * (e * i_elem - f * h) - b * (d * i_elem - f * g) + c * (d * h - e * g);
                    if det <= 0.0 {
                        return Err(SklearsError::NumericalError(format!(
                            "Matrix {} has non-positive determinant",
                            i
                        )));
                    }
                    log_dets[i] = det.ln();
                }
            }
            _ => {
                // General case - would use LU decomposition
                return Err(SklearsError::InvalidInput(
                    "Only 2x2 and 3x3 matrices supported in batch mode".to_string(),
                ));
            }
        }

        Ok(log_dets)
    }

    /// SIMD-accelerated quadratic form computation: x^T A x for multiple vectors
    pub fn simd_batch_quadratic_form(
        &self,
        vectors: &Array2<Float>,
        matrix: &Array2<Float>,
    ) -> Result<Array1<Float>> {
        let n_vectors = vectors.nrows();
        let mut results = Array1::zeros(n_vectors);

        if self.config.parallel_simd && n_vectors >= 64 {
            let parallel_results: Result<Vec<Float>> = (0..n_vectors)
                .into_par_iter()
                .map(|i| -> Result<Float> {
                    let x = vectors.row(i);
                    let temp = self.base_ops.simd_matvec(matrix, &x.to_owned())?;
                    self.base_ops.simd_dot_product(&x, temp.view())
                })
                .collect();

            let computed_results = parallel_results?;
            for (i, value) in computed_results.into_iter().enumerate() {
                results[i] = value;
            }
        } else {
            for i in 0..n_vectors {
                let x = vectors.row(i);
                let temp = self.base_ops.simd_matvec(matrix, &x.to_owned())?;
                results[i] = self.base_ops.simd_dot_product(&x, temp.view())?;
            }
        }

        Ok(results)
    }

    /// SIMD-accelerated distance matrix computation
    pub fn simd_pairwise_distances(
        &self,
        x: &Array2<Float>,
        y: &Array2<Float>,
    ) -> Result<Array2<Float>> {
        let (n_x, dim_x) = x.dim();
        let (n_y, dim_y) = y.dim();

        if dim_x != dim_y {
            return Err(SklearsError::InvalidInput(
                "X and Y must have same number of features".to_string(),
            ));
        }

        let mut distances = Array2::zeros((n_x, n_y));

        if self.config.parallel_simd {
            let results: Result<Vec<Vec<Float>>> = (0..n_x)
                .into_par_iter()
                .map(|i| -> Result<Vec<Float>> {
                    let x_row = x.row(i);
                    let mut row_distances = Vec::with_capacity(n_y);

                    for j in 0..n_y {
                        let y_row = y.row(j);
                        let diff = self
                            .base_ops
                            .simd_element_wise_subtract(&x_row.to_owned(), &y_row.to_owned())?;
                        let dist_sq = self.base_ops.simd_dot_product(&diff.view(), diff.view())?;
                        row_distances.push(dist_sq.sqrt());
                    }
                    Ok(row_distances)
                })
                .collect();

            let computed_results = results?;
            for (i, row_data) in computed_results.into_iter().enumerate() {
                for (j, dist) in row_data.into_iter().enumerate() {
                    distances[[i, j]] = dist;
                }
            }
        } else {
            for i in 0..n_x {
                for j in 0..n_y {
                    let x_row = x.row(i);
                    let y_row = y.row(j);
                    let diff = self
                        .base_ops
                        .simd_element_wise_subtract(&x_row.to_owned(), &y_row.to_owned())?;
                    let dist_sq = self.base_ops.simd_dot_product(&diff.view(), diff.view())?;
                    distances[[i, j]] = dist_sq.sqrt();
                }
            }
        }

        Ok(distances)
    }

    /// SIMD-accelerated cross-covariance computation
    pub fn simd_cross_covariance(
        &self,
        x: &Array2<Float>,
        y: &Array2<Float>,
    ) -> Result<Array2<Float>> {
        if x.nrows() != y.nrows() {
            return Err(SklearsError::InvalidInput(
                "X and Y must have same number of samples".to_string(),
            ));
        }

        let n_samples = x.nrows() as Float;

        // Compute means
        let mean_x = self.base_ops.simd_column_mean(x)?;
        let mean_y = self.base_ops.simd_column_mean(y)?;

        // Center the data
        let mut centered_x = Array2::zeros(x.raw_dim());
        let mut centered_y = Array2::zeros(y.raw_dim());

        for (i, (mut cx_row, mut cy_row)) in centered_x
            .axis_iter_mut(Axis(0))
            .zip(centered_y.axis_iter_mut(Axis(0)))
            .enumerate()
        {
            let x_row = x.row(i);
            let y_row = y.row(i);
            cx_row.assign(
                &self
                    .base_ops
                    .simd_element_wise_subtract(&x_row.to_owned(), &mean_x)?,
            );
            cy_row.assign(
                &self
                    .base_ops
                    .simd_element_wise_subtract(&y_row.to_owned(), &mean_y)?,
            );
        }

        // Compute cross-covariance: X^T * Y / (n - 1)
        let cross_cov = self
            .base_ops
            .simd_matmul(&centered_x.t().to_owned(), &centered_y)?;
        Ok(cross_cov / (n_samples - 1.0))
    }

    /// SIMD-accelerated feature scaling (standardization)
    pub fn simd_standardize(&self, data: &Array2<Float>) -> Result<Array2<Float>> {
        let mean = self.base_ops.simd_column_mean(data)?;

        // Compute standard deviation for each feature
        let mut std_dev = Array1::zeros(data.ncols());
        for (j, std_elem) in std_dev.iter_mut().enumerate() {
            let col = data.column(j);
            let mean_j = mean[j];
            let variance = col.iter().map(|&x| (x - mean_j).powi(2)).sum::<Float>()
                / (data.nrows() - 1) as Float;
            *std_elem = variance.sqrt();
        }

        // Standardize: (x - mean) / std
        let mut standardized = Array2::zeros(data.raw_dim());
        for (i, mut std_row) in standardized.axis_iter_mut(Axis(0)).enumerate() {
            let data_row = data.row(i);
            let centered = self
                .base_ops
                .simd_element_wise_subtract(&data_row.to_owned(), &mean)?;
            let scaled = self.element_wise_divide(&centered, &std_dev)?;
            std_row.assign(&scaled);
        }

        Ok(standardized)
    }

    /// Element-wise division helper
    fn element_wise_divide(&self, a: &Array1<Float>, b: &Array1<Float>) -> Result<Array1<Float>> {
        if a.len() != b.len() {
            return Err(SklearsError::InvalidInput(
                "Arrays must have same length".to_string(),
            ));
        }

        Ok(a.iter()
            .zip(b.iter())
            .map(|(&a_val, &b_val)| {
                if b_val.abs() < Float::EPSILON {
                    0.0 // Handle division by zero
                } else {
                    a_val / b_val
                }
            })
            .collect::<Array1<Float>>())
    }

    /// SIMD-accelerated log-sum-exp computation for numerical stability
    pub fn simd_log_sum_exp(&self, x: &Array1<Float>) -> Result<Float> {
        if x.is_empty() {
            return Err(SklearsError::InvalidInput(
                "Input array is empty".to_string(),
            ));
        }

        // Find maximum for numerical stability
        let max_val = x.iter().fold(Float::NEG_INFINITY, |acc, &val| acc.max(val));

        // Compute log(sum(exp(x - max))) + max
        let shifted = x.mapv(|val| val - max_val);
        let exp_values = self.simd_exp(&shifted)?;
        let sum_exp = exp_values.sum();

        Ok(sum_exp.ln() + max_val)
    }

    /// Performance profiling for SIMD operations
    pub fn benchmark_simd_vs_scalar(&self, size: usize) -> SimdBenchmarkResults {
        let a = Array1::from_vec((0..size).map(|i| i as Float).collect());
        let b = Array1::from_vec((0..size).map(|i| (i + 1) as Float).collect());

        // SIMD benchmark
        let start = std::time::Instant::now();
        let _simd_result = self.base_ops.simd_dot_product(&a.view(), b.view());
        let simd_time = start.elapsed();

        // Scalar benchmark
        let start = std::time::Instant::now();
        let _scalar_result: Float = a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum();
        let scalar_time = start.elapsed();

        SimdBenchmarkResults {
            size,
            simd_time_ns: simd_time.as_nanos() as u64,
            scalar_time_ns: scalar_time.as_nanos() as u64,
            speedup: scalar_time.as_nanos() as f64 / simd_time.as_nanos() as f64,
        }
    }
}

impl Default for AdvancedSimdOps {
    fn default() -> Self {
        Self::new()
    }
}

/// Benchmark results comparing SIMD vs scalar performance
#[derive(Debug, Clone)]
pub struct SimdBenchmarkResults {
    /// size
    pub size: usize,
    /// simd_time_ns
    pub simd_time_ns: u64,
    /// scalar_time_ns
    pub scalar_time_ns: u64,
    /// speedup
    pub speedup: f64,
}

/// ARM NEON SIMD support (placeholder for future ARM optimization)
#[cfg(target_arch = "aarch64")]
pub struct NeonSimdOps {
    config: SimdConfig,
}

#[cfg(target_arch = "aarch64")]
impl Default for NeonSimdOps {
    fn default() -> Self {
        Self::new()
    }
}

#[cfg(target_arch = "aarch64")]
impl NeonSimdOps {
    pub fn new() -> Self {
        Self {
            config: SimdConfig::default(),
        }
    }

    /// ARM NEON dot product implementation (placeholder)
    pub fn neon_dot_product(&self, a: &[Float], b: &[Float]) -> Result<Float> {
        // Placeholder for ARM NEON implementation
        // Real implementation would use ARM NEON intrinsics
        Ok(a.iter().zip(b.iter()).map(|(&x, &y)| x * y).sum())
    }

    /// Check NEON support
    pub fn check_neon_support(&self) -> bool {
        // Placeholder for ARM NEON detection
        true
    }
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;
    use scirs2_core::ndarray::array;

    #[test]
    fn test_simd_dot_product() {
        let simd_ops = SimdMatrixOps::new();
        let a = array![1.0, 2.0, 3.0, 4.0];
        let b = array![5.0, 6.0, 7.0, 8.0];

        let result = simd_ops.simd_dot_product(&a.view(), b.view()).unwrap();
        let expected = 1.0 * 5.0 + 2.0 * 6.0 + 3.0 * 7.0 + 4.0 * 8.0; // = 70.0

        assert_abs_diff_eq!(result, expected, epsilon = 1e-10);
    }

    #[test]
    fn test_simd_matvec() {
        let simd_ops = SimdMatrixOps::new();
        let matrix = array![[1.0, 2.0], [3.0, 4.0]];
        let vector = array![5.0, 6.0];

        let result = simd_ops.simd_matvec(&matrix, &vector).unwrap();
        let expected = array![1.0 * 5.0 + 2.0 * 6.0, 3.0 * 5.0 + 4.0 * 6.0]; // = [17.0, 39.0]

        assert_abs_diff_eq!(result[0], expected[0], epsilon = 1e-10);
        assert_abs_diff_eq!(result[1], expected[1], epsilon = 1e-10);
    }

    #[test]
    fn test_simd_matmul() {
        let simd_ops = SimdMatrixOps::new();
        let a = array![[1.0, 2.0], [3.0, 4.0]];
        let b = array![[5.0, 6.0], [7.0, 8.0]];

        let result = simd_ops.simd_matmul(&a, &b).unwrap();
        // Expected: [[19.0, 22.0], [43.0, 50.0]]

        assert_abs_diff_eq!(result[[0, 0]], 19.0, epsilon = 1e-10);
        assert_abs_diff_eq!(result[[0, 1]], 22.0, epsilon = 1e-10);
        assert_abs_diff_eq!(result[[1, 0]], 43.0, epsilon = 1e-10);
        assert_abs_diff_eq!(result[[1, 1]], 50.0, epsilon = 1e-10);
    }

    #[test]
    fn test_simd_element_wise_ops() {
        let simd_ops = SimdMatrixOps::new();
        let a = array![1.0, 2.0, 3.0, 4.0];
        let b = array![5.0, 6.0, 7.0, 8.0];

        let add_result = simd_ops.simd_element_wise_add(&a, &b).unwrap();
        let expected_add = array![6.0, 8.0, 10.0, 12.0];

        for (r, e) in add_result.iter().zip(expected_add.iter()) {
            assert_abs_diff_eq!(*r, *e, epsilon = 1e-10);
        }

        let sub_result = simd_ops.simd_element_wise_subtract(&a, &b).unwrap();
        let expected_sub = array![-4.0, -4.0, -4.0, -4.0];

        for (r, e) in sub_result.iter().zip(expected_sub.iter()) {
            assert_abs_diff_eq!(*r, *e, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_simd_covariance() {
        let simd_ops = SimdMatrixOps::new();
        let data = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]];

        let cov = simd_ops.simd_covariance(&data).unwrap();

        // Expected covariance matrix for this simple case
        assert!(cov.nrows() == 2 && cov.ncols() == 2);
        assert!(cov[[0, 0]] > 0.0); // Variance should be positive
        assert!(cov[[1, 1]] > 0.0); // Variance should be positive
        assert_abs_diff_eq!(cov[[0, 1]], cov[[1, 0]], epsilon = 1e-10); // Should be symmetric
    }

    #[test]
    fn test_simd_support_detection() {
        let simd_ops = SimdMatrixOps::new();
        let support = simd_ops.check_simd_support();

        // Just verify the function runs without panicking
        println!("SIMD Support: {:?}", support);
    }

    #[test]
    fn test_simd_array_ops_trait() {
        let a = array![1.0, 2.0, 3.0, 4.0];
        let b = array![5.0, 6.0, 7.0, 8.0];

        let dot_result = a.simd_dot(&b).unwrap();
        let expected = 70.0;
        assert_abs_diff_eq!(dot_result, expected, epsilon = 1e-10);

        let add_result = a.simd_add(&b).unwrap();
        let expected_add = array![6.0, 8.0, 10.0, 12.0];

        for (r, e) in add_result.iter().zip(expected_add.iter()) {
            assert_abs_diff_eq!(*r, *e, epsilon = 1e-10);
        }
    }
}