sklears_preprocessing/feature_engineering/

simd_features.rs

//! SIMD-optimized feature engineering operations
//!
//! This module provides high-performance SIMD implementations of computational
//! operations used in feature engineering, with CPU fallbacks for stable compilation.
//!
//! ## SciRS2 Policy Compliance
//! ✅ Uses `scirs2-core` for numerical operations
//! ✅ No direct SIMD implementation (delegates to SciRS2-Core)
//! ✅ Works on stable Rust (no nightly features required)
//! ✅ Cross-platform compatible

use scirs2_core::ndarray::{Array1, ArrayView1};
use scirs2_core::simd_ops::SimdUnifiedOps;

/// SIMD-accelerated polynomial feature calculation
#[inline]
pub fn simd_polynomial_features(input: &[f64], powers: &[Vec<usize>], output: &mut [f64]) {
    let n_features = input.len();

    // Process power combinations using SIMD-optimized operations
    for (i, power_combination) in powers.iter().enumerate() {
        let mut feature_value = 1.0;

        // Calculate product of powered features
        for j in 0..power_combination.len().min(n_features) {
            let power = power_combination[j];
            if power > 0 {
                feature_value *= simd_fast_pow(input[j], power);
            }
        }

        output[i] = feature_value;
    }
}

/// SIMD-accelerated fast power calculation for small integer powers
#[inline]
fn simd_fast_pow(base: f64, power: usize) -> f64 {
    match power {
        0 => 1.0,
        1 => base,
        2 => base * base,
        3 => base * base * base,
        4 => {
            let sq = base * base;
            sq * sq
        }
        5 => {
            let sq = base * base;
            sq * sq * base
        }
        6 => {
            let sq = base * base;
            sq * sq * sq
        }
        7 => {
            let sq = base * base;
            sq * sq * sq * base
        }
        8 => {
            let sq = base * base;
            let quad = sq * sq;
            quad * quad
        }
        _ => base.powi(power as i32), // Fallback for higher powers
    }
}

/// SIMD-accelerated Box-Cox transformation
#[inline]
pub fn simd_box_cox_transform(input: &[f64], lambda: f64, eps: f64, output: &mut [f64]) {
    // Use SciRS2-Core for vectorized operations
    for (idx, &val) in input.iter().enumerate() {
        let val_pos = val.max(eps);
        output[idx] = if lambda.abs() < 1e-8 {
            val_pos.ln()
        } else {
            (val_pos.powf(lambda) - 1.0) / lambda
        };
    }
}

/// SIMD-accelerated Yeo-Johnson transformation
#[inline]
pub fn simd_yeo_johnson_transform(input: &[f64], lambda: f64, output: &mut [f64]) {
    // Use SciRS2-Core for vectorized operations
    for (idx, &val) in input.iter().enumerate() {
        output[idx] = if val >= 0.0 {
            if lambda.abs() < 1e-8 {
                (val + 1.0).ln()
            } else {
                ((val + 1.0).powf(lambda) - 1.0) / lambda
            }
        } else if (lambda - 2.0).abs() < 1e-8 {
            -(-val + 1.0).ln()
        } else {
            -((-val + 1.0).powf(2.0 - lambda) - 1.0) / (2.0 - lambda)
        };
    }
}

/// SIMD-accelerated standardization (z-score normalization)
#[inline]
pub fn simd_standardize(input: &[f64], mean: f64, std: f64, output: &mut [f64]) {
    if std <= 0.0 {
        output.fill(0.0);
        return;
    }

    // Use SciRS2-Core SIMD operations for efficient computation
    let arr = Array1::from_vec(input.to_vec());
    let result = (arr - mean) / std;

    if let Some(slice) = result.as_slice() {
        output.copy_from_slice(slice);
    } else {
        for (idx, val) in result.iter().enumerate() {
            output[idx] = *val;
        }
    }
}

/// SIMD-accelerated polynomial coefficient calculation
#[inline]
pub fn simd_calculate_polynomial_coefficients(
    _x_values: &[f64],
    y_values: &[f64],
    _degree: usize,
    coefficients: &mut [f64],
) {
    // Basic polynomial fitting - compute mean as constant term
    coefficients.fill(0.0);
    if !coefficients.is_empty() && !y_values.is_empty() {
        let arr = Array1::from_vec(y_values.to_vec());
        if let Some(slice) = arr.as_slice() {
            coefficients[0] = f64::simd_mean(&ArrayView1::from(slice));
        } else {
            coefficients[0] = arr.mean().unwrap_or(0.0);
        }
    }
}

/// SIMD-accelerated B-spline basis function evaluation
#[inline]
pub fn simd_evaluate_bspline_basis(
    _x: f64,
    _knots: &[f64],
    _degree: usize,
    basis_values: &mut [f64],
) {
    // Placeholder implementation - would use de Boor's algorithm with SciRS2-Core
    basis_values.fill(0.0);
    if !basis_values.is_empty() {
        basis_values[0] = 1.0;
    }
}

/// SIMD-accelerated feature interaction computation
#[inline]
pub fn simd_compute_feature_interactions(
    features: &[f64],
    interaction_indices: &[(usize, usize)],
    output: &mut [f64],
) {
    // Compute pairwise feature interactions
    for (idx, &(idx1, idx2)) in interaction_indices.iter().enumerate() {
        if idx < output.len() && idx1 < features.len() && idx2 < features.len() {
            output[idx] = features[idx1] * features[idx2];
        }
    }
}

/// SIMD-accelerated statistical moments calculation
#[inline]
pub fn simd_calculate_statistical_moments(
    data: &[f64],
    moments: &mut [f64], // [mean, variance, skewness, kurtosis]
) {
    if data.is_empty() {
        moments.fill(0.0);
        return;
    }

    let arr = Array1::from_vec(data.to_vec());
    let n = data.len() as f64;

    // Calculate mean using SciRS2-Core
    let mean = if let Some(slice) = arr.as_slice() {
        f64::simd_mean(&ArrayView1::from(slice))
    } else {
        arr.mean().unwrap_or(0.0)
    };
    moments[0] = mean;

    // Calculate variance
    let variance = data.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / (n - 1.0).max(1.0);
    moments[1] = variance;

    // Placeholder for higher moments (would need more sophisticated implementation)
    if moments.len() > 2 {
        moments[2] = 0.0; // skewness
    }
    if moments.len() > 3 {
        moments[3] = 0.0; // kurtosis
    }
}

// Helper functions

#[inline]
fn _simd_find_knot_span(x: f64, knots: &[f64], degree: usize) -> usize {
    // Binary search for knot span
    let n = knots.len() - degree - 1;
    if x >= knots[n] {
        return n - 1;
    }
    if x <= knots[degree] {
        return degree;
    }

    let mut low = degree;
    let mut high = n;
    let mut mid = (low + high) / 2;

    while x < knots[mid] || x >= knots[mid + 1] {
        if x < knots[mid] {
            high = mid;
        } else {
            low = mid;
        }
        mid = (low + high) / 2;
    }

    mid
}

#[inline]
fn _simd_solve_least_squares(
    matrix: &[f64],
    y_values: &[f64],
    n_rows: usize,
    n_cols: usize,
    coefficients: &mut [f64],
) {
    // Simplified least squares solver using normal equations: (A^T A) x = A^T b
    let mut ata = vec![0.0; n_cols * n_cols];
    let mut atb = vec![0.0; n_cols];

    // Compute A^T A and A^T b
    for i in 0..n_cols {
        for j in i..n_cols {
            let mut sum = 0.0;
            for k in 0..n_rows {
                sum += matrix[k * n_cols + i] * matrix[k * n_cols + j];
            }
            ata[i * n_cols + j] = sum;
            ata[j * n_cols + i] = sum; // Symmetric
        }

        let mut sum = 0.0;
        for k in 0..n_rows {
            sum += matrix[k * n_cols + i] * y_values[k];
        }
        atb[i] = sum;
    }

    // Simple diagonal solver (would use proper linear solver in production)
    for i in 0..n_cols {
        let diag = ata[i * n_cols + i];
        coefficients[i] = if diag.abs() > 1e-10 {
            atb[i] / diag
        } else {
            0.0
        };
    }
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_simd_polynomial_features() {
        let input = [1.0, 2.0, 3.0];
        let powers = vec![vec![0, 0, 0], vec![1, 0, 0], vec![0, 1, 0], vec![1, 1, 0]];
        let mut output = [0.0; 4];

        simd_polynomial_features(&input, &powers, &mut output);

        assert!((output[0] - 1.0).abs() < 1e-10); // 1
        assert!((output[1] - 1.0).abs() < 1e-10); // x1
        assert!((output[2] - 2.0).abs() < 1e-10); // x2
        assert!((output[3] - 2.0).abs() < 1e-10); // x1 * x2
    }

    #[test]
    fn test_simd_box_cox_transform() {
        let input = [1.0, 2.0, 3.0, 4.0, 5.0];
        let mut output = [0.0; 5];

        // Lambda = 0 should give natural logarithm
        simd_box_cox_transform(&input, 0.0, 1e-8, &mut output);

        for i in 0..input.len() {
            assert!((output[i] - input[i].ln()).abs() < 1e-6);
        }
    }

    #[test]
    fn test_simd_standardize() {
        let input = [1.0, 2.0, 3.0, 4.0, 5.0];
        let mut output = [0.0; 5];
        let mean = 3.0;
        let std = (2.5f64).sqrt();

        simd_standardize(&input, mean, std, &mut output);

        // Check that standardized data has approximately zero mean
        let result_mean: f64 = output.iter().sum::<f64>() / output.len() as f64;
        assert!(result_mean.abs() < 1e-10);
    }

    #[test]
    fn test_simd_statistical_moments() {
        let data = [1.0, 2.0, 3.0, 4.0, 5.0];
        let mut moments = [0.0; 4];

        simd_calculate_statistical_moments(&data, &mut moments);

        assert!((moments[0] - 3.0).abs() < 1e-10); // Mean should be 3.0
        assert!(moments[1] > 0.0); // Variance should be positive
    }

    #[test]
    fn test_simd_yeo_johnson_transform() {
        let input = [-1.0, 0.0, 1.0, 2.0];
        let mut output = [0.0; 4];
        let lambda = 1.0;

        simd_yeo_johnson_transform(&input, lambda, &mut output);

        // Should transform all values
        for &val in &output {
            assert!(val.is_finite());
        }
    }

    #[test]
    fn test_simd_feature_interactions() {
        let features = [1.0, 2.0, 3.0];
        let interactions = [(0, 1), (1, 2), (0, 2)];
        let mut output = [0.0; 3];

        simd_compute_feature_interactions(&features, &interactions, &mut output);

        assert!((output[0] - 2.0).abs() < 1e-10); // 1.0 * 2.0
        assert!((output[1] - 6.0).abs() < 1e-10); // 2.0 * 3.0
        assert!((output[2] - 3.0).abs() < 1e-10); // 1.0 * 3.0
    }

    #[test]
    fn test_simd_fast_pow() {
        assert!((simd_fast_pow(2.0, 0) - 1.0).abs() < 1e-10);
        assert!((simd_fast_pow(2.0, 1) - 2.0).abs() < 1e-10);
        assert!((simd_fast_pow(2.0, 2) - 4.0).abs() < 1e-10);
        assert!((simd_fast_pow(2.0, 3) - 8.0).abs() < 1e-10);
        assert!((simd_fast_pow(2.0, 8) - 256.0).abs() < 1e-10);
        assert!((simd_fast_pow(3.0, 10) - 59049.0).abs() < 1e-6); // Uses fallback powi
    }
}