rusty_backtest/

simd_enhanced.rs

//! Enhanced SIMD Operations for High-Frequency Trading
//!
//! This module provides advanced SIMD optimizations specifically designed for HFT workloads:
//! - Cache-aligned memory layouts
//! - Optimized memory access patterns
//! - Branch-free algorithms where possible
//! - Prefetching for large datasets
//! - Structure-of-Arrays (SoA) optimizations

#[cfg(test)]
use std::hint::black_box;
use wide::f64x4;

/// Portable prefetch function that compiles to optimal instructions per platform
#[inline(always)]
fn prefetch_data(ptr: *const f64) {
    #[cfg(target_arch = "x86_64")]
    {
        unsafe {
            std::arch::x86_64::_mm_prefetch(ptr as *const i8, std::arch::x86_64::_MM_HINT_T0);
        }
    }
    #[cfg(target_arch = "aarch64")]
    {
        unsafe {
            std::arch::aarch64::_prefetch(
                ptr as *const i8,
                std::arch::aarch64::_PREFETCH_READ,
                std::arch::aarch64::_PREFETCH_LOCALITY3,
            );
        }
    }
    // On other architectures, prefetch is a no-op
    // The compiler may still optimize memory access patterns
    #[cfg(not(any(target_arch = "x86_64", target_arch = "aarch64")))]
    {
        let _ = ptr; // Suppress unused variable warning
    }
}

/// Re-export the unified 64-byte aligned buffer type
///
/// This type alias maintains backward compatibility while using the
/// consolidated aligned buffer implementation from rusty-common.
pub use rusty_common::AlignedBuffer64 as SimdAlignedBuffer;

/// Enhanced SIMD operations optimized for financial data processing
pub struct EnhancedSimdOps;

impl EnhancedSimdOps {
    /// Optimized SIMD sum with explicit prefetching for large arrays
    #[inline]
    pub fn sum_f64_optimized(values: &[f64]) -> f64 {
        if values.is_empty() {
            return 0.0;
        }

        if values.len() == 1 {
            return values[0];
        }

        if values.len() < 4 {
            return values.iter().sum();
        }

        // For large arrays, use prefetching
        if values.len() > 1024 {
            Self::sum_f64_with_prefetch(values)
        } else {
            Self::sum_f64_fast(values)
        }
    }

    /// Fast SIMD sum for medium-sized arrays without prefetching overhead
    #[inline]
    fn sum_f64_fast(values: &[f64]) -> f64 {
        let mut i = 0;
        let len = values.len();

        // Use two accumulators to reduce dependency chains
        let mut sum1 = f64x4::ZERO;
        let mut sum2 = f64x4::ZERO;

        // Process 8 elements at a time (2 SIMD lanes)
        while i + 8 <= len {
            let v1 = f64x4::from([values[i], values[i + 1], values[i + 2], values[i + 3]]);
            let v2 = f64x4::from([values[i + 4], values[i + 5], values[i + 6], values[i + 7]]);

            sum1 += v1;
            sum2 += v2;

            i += 8;
        }

        // Process remaining 4-element chunks
        while i + 4 <= len {
            let v = f64x4::from([values[i], values[i + 1], values[i + 2], values[i + 3]]);
            sum1 += v;
            i += 4;
        }

        // Combine the two accumulators
        let combined = sum1 + sum2;
        let combined_array = combined.as_array_ref();
        let mut total =
            combined_array[0] + combined_array[1] + combined_array[2] + combined_array[3];

        // Handle remaining elements
        while i < len {
            total += values[i];
            i += 1;
        }

        total
    }

    /// SIMD sum with explicit prefetching for large arrays
    #[inline]
    fn sum_f64_with_prefetch(values: &[f64]) -> f64 {
        let mut i = 0;
        let len = values.len();

        // Use four accumulators to maximize instruction-level parallelism
        let mut sum1 = f64x4::ZERO;
        let mut sum2 = f64x4::ZERO;
        let mut sum3 = f64x4::ZERO;
        let mut sum4 = f64x4::ZERO;

        // Process 16 elements at a time with prefetching
        while i + 16 <= len {
            // Prefetch next cache line (64 bytes = 8 f64 values ahead)
            if i + 24 < len {
                unsafe {
                    prefetch_data(values.as_ptr().add(i + 16));
                }
            }

            let v1 = f64x4::from([values[i], values[i + 1], values[i + 2], values[i + 3]]);
            let v2 = f64x4::from([values[i + 4], values[i + 5], values[i + 6], values[i + 7]]);
            let v3 = f64x4::from([values[i + 8], values[i + 9], values[i + 10], values[i + 11]]);
            let v4 = f64x4::from([
                values[i + 12],
                values[i + 13],
                values[i + 14],
                values[i + 15],
            ]);

            sum1 += v1;
            sum2 += v2;
            sum3 += v3;
            sum4 += v4;

            i += 16;
        }

        // Combine all accumulators
        let combined = (sum1 + sum2) + (sum3 + sum4);
        let combined_array = combined.as_array_ref();
        let mut total =
            combined_array[0] + combined_array[1] + combined_array[2] + combined_array[3];

        // Handle remaining elements with scalar code
        while i < len {
            total += values[i];
            i += 1;
        }

        total
    }

    /// Welford's algorithm for variance calculation (scalar implementation)
    /// This is the actual Welford's algorithm with single-pass, numerically stable computation
    #[inline]
    pub fn variance_welford_scalar(values: &[f64]) -> f64 {
        if values.len() < 2 {
            return 0.0;
        }

        // Initialize for Welford's algorithm
        let mut count = 0f64;
        let mut mean = 0.0;
        let mut m2 = 0.0;

        // Process elements using Welford's update formula
        for &value in values {
            count += 1.0;
            let delta = value - mean;
            mean += delta / count;
            let delta2 = value - mean;
            m2 += delta * delta2;
        }

        m2 / count
    }

    /// Two-pass variance calculation using SIMD with aligned memory access
    /// Uses simd_aligned for guaranteed memory alignment and optimal performance
    #[inline]
    pub fn variance_two_pass_simd(values: &[f64]) -> f64 {
        if values.len() < 2 {
            return 0.0;
        }

        let len = values.len();
        let len_f64 = len as f64;

        // Use simple scalar calculation for small arrays where SIMD overhead isn't worth it
        if len <= 4 {
            return Self::variance_welford_scalar(values);
        }

        // Calculate mean using SIMD on aligned buffer but only for complete chunks
        let complete_chunks = len / 4;
        let mut sum = f64x4::splat(0.0);

        // Process complete 4-element chunks using SIMD
        for i in 0..complete_chunks {
            let start_idx = i * 4;
            let chunk = f64x4::from([
                values[start_idx],
                values[start_idx + 1],
                values[start_idx + 2],
                values[start_idx + 3],
            ]);
            sum += chunk;
        }

        // Sum the SIMD lanes
        let sum_array = sum.as_array_ref();
        let mut total_sum = sum_array[0] + sum_array[1] + sum_array[2] + sum_array[3];

        // Add remaining elements using scalar arithmetic
        for i in (complete_chunks * 4)..len {
            total_sum += values[i];
        }

        let mean = total_sum / len_f64;

        // Second pass: Calculate variance using SIMD for complete chunks
        let mean_vec = f64x4::splat(mean);
        let mut variance_sum = f64x4::splat(0.0);

        // Process complete 4-element chunks using SIMD
        for i in 0..complete_chunks {
            let start_idx = i * 4;
            let chunk = f64x4::from([
                values[start_idx],
                values[start_idx + 1],
                values[start_idx + 2],
                values[start_idx + 3],
            ]);
            let diff = chunk - mean_vec;
            variance_sum += diff * diff;
        }

        // Sum the SIMD lanes
        let variance_array = variance_sum.as_array_ref();
        let mut total_variance =
            variance_array[0] + variance_array[1] + variance_array[2] + variance_array[3];

        // Add remaining elements using scalar arithmetic
        for i in (complete_chunks * 4)..len {
            let diff = values[i] - mean;
            total_variance += diff * diff;
        }

        total_variance / len_f64
    }

    /// Calculate running averages for price series using SIMD
    /// Optimized for typical HFT window sizes (50-200 periods)
    #[inline]
    pub fn running_average_simd(values: &[f64], window_size: usize) -> Vec<f64> {
        if values.len() < window_size {
            return vec![];
        }

        let mut results = Vec::with_capacity(values.len() - window_size + 1);

        // Calculate first window sum
        let mut window_sum = Self::sum_f64_fast(&values[0..window_size]);
        results.push(window_sum / window_size as f64);

        // Use sliding window technique for subsequent calculations
        for i in window_size..values.len() {
            window_sum = window_sum - values[i - window_size] + values[i];
            results.push(window_sum / window_size as f64);
        }

        results
    }

    /// Optimized correlation calculation using SIMD
    /// Processes two price series to calculate correlation coefficient
    #[inline]
    pub fn correlation_simd(x: &[f64], y: &[f64]) -> f64 {
        if x.len() != y.len() || x.len() < 2 {
            return 0.0;
        }

        let n = x.len() as f64;
        let sum_x = Self::sum_f64_optimized(x);
        let sum_y = Self::sum_f64_optimized(y);
        let mean_x = sum_x / n;
        let mean_y = sum_y / n;

        let mut sum_xy = 0.0;
        let mut sum_x2 = 0.0;
        let mut sum_y2 = 0.0;

        let mut i = 0;
        let len = x.len();

        // Process 4 elements at a time
        while i + 4 <= len {
            let x_vec = f64x4::from([x[i], x[i + 1], x[i + 2], x[i + 3]]);
            let y_vec = f64x4::from([y[i], y[i + 1], y[i + 2], y[i + 3]]);
            let mean_x_vec = f64x4::splat(mean_x);
            let mean_y_vec = f64x4::splat(mean_y);

            let dx = x_vec - mean_x_vec;
            let dy = y_vec - mean_y_vec;

            let xy = dx * dy;
            let x2 = dx * dx;
            let y2 = dy * dy;

            let xy_array = xy.as_array_ref();
            let x2_array = x2.as_array_ref();
            let y2_array = y2.as_array_ref();

            sum_xy += xy_array[0] + xy_array[1] + xy_array[2] + xy_array[3];
            sum_x2 += x2_array[0] + x2_array[1] + x2_array[2] + x2_array[3];
            sum_y2 += y2_array[0] + y2_array[1] + y2_array[2] + y2_array[3];

            i += 4;
        }

        // Handle remaining elements
        while i < len {
            let dx = x[i] - mean_x;
            let dy = y[i] - mean_y;
            sum_xy += dx * dy;
            sum_x2 += dx * dx;
            sum_y2 += dy * dy;
            i += 1;
        }

        let denominator = (sum_x2 * sum_y2).sqrt();
        if denominator > f64::EPSILON {
            sum_xy / denominator
        } else {
            0.0
        }
    }

    /// Optimized price returns calculation with SIMD
    /// Calculates log returns: ln(price\[i\] / price\[i-1\])
    #[inline]
    pub fn log_returns_simd(prices: &[f64]) -> Vec<f64> {
        if prices.len() < 2 {
            return vec![];
        }

        let mut returns = Vec::with_capacity(prices.len() - 1);

        let mut i = 1;
        let len = prices.len();

        // Process 4 returns at a time
        while i + 4 <= len {
            let curr_prices = f64x4::from([prices[i], prices[i + 1], prices[i + 2], prices[i + 3]]);
            let prev_prices = f64x4::from([prices[i - 1], prices[i], prices[i + 1], prices[i + 2]]);

            // Calculate ratios and take natural log
            let ratios = curr_prices / prev_prices;
            let ratios_array = ratios.as_array_ref();

            returns.push(ratios_array[0].ln());
            returns.push(ratios_array[1].ln());
            returns.push(ratios_array[2].ln());
            returns.push(ratios_array[3].ln());

            i += 4;
        }

        // Handle remaining elements
        while i < len {
            let log_return = (prices[i] / prices[i - 1]).ln();
            returns.push(log_return);
            i += 1;
        }

        returns
    }

    /// Structure-of-Arrays optimized order book level processing
    /// Processes bid/ask levels more efficiently by separating prices and quantities
    #[inline]
    pub fn process_order_book_levels_soa(
        prices: &[f64],
        quantities: &[f64],
        operation: impl Fn(f64, f64) -> f64,
    ) -> Vec<f64> {
        if prices.len() != quantities.len() {
            return vec![];
        }

        let mut results = Vec::with_capacity(prices.len());
        let mut i = 0;
        let len = prices.len();

        // Process 4 levels at a time
        while i + 4 <= len {
            let price_vec = f64x4::from([prices[i], prices[i + 1], prices[i + 2], prices[i + 3]]);
            let qty_vec = f64x4::from([
                quantities[i],
                quantities[i + 1],
                quantities[i + 2],
                quantities[i + 3],
            ]);

            let price_array = price_vec.as_array_ref();
            let qty_array = qty_vec.as_array_ref();

            for j in 0..4 {
                #[cfg(test)]
                let result = black_box(operation(price_array[j], qty_array[j]));
                #[cfg(not(test))]
                let result = operation(price_array[j], qty_array[j]);
                results.push(result);
            }

            i += 4;
        }

        // Handle remaining elements
        while i < len {
            results.push(operation(prices[i], quantities[i]));
            i += 1;
        }

        results
    }
}

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

    #[test]
    fn test_enhanced_sum_small_array() {
        let values = [1.0, 2.0, 3.0];
        let sum = EnhancedSimdOps::sum_f64_optimized(&values);
        assert_eq!(sum, 6.0);
    }

    #[test]
    fn test_enhanced_sum_medium_array() {
        let values: Vec<f64> = (1..=100).map(|i| i as f64).collect();
        let sum = EnhancedSimdOps::sum_f64_optimized(&values);
        assert_eq!(sum, 5050.0);
    }

    #[test]
    fn test_enhanced_sum_large_array() {
        let values: Vec<f64> = (1..=10000).map(|i| i as f64).collect();
        let sum = EnhancedSimdOps::sum_f64_optimized(&values);
        assert_eq!(sum, 50005000.0);
    }

    #[test]
    fn test_variance_welford() {
        let values = [2.0, 4.0, 6.0, 8.0];
        let variance = EnhancedSimdOps::variance_welford_scalar(&values);
        assert!((variance - 5.0).abs() < 1e-10);
    }

    #[test]
    fn test_running_average() {
        let values = [1.0, 2.0, 3.0, 4.0, 5.0];
        let averages = EnhancedSimdOps::running_average_simd(&values, 3);
        assert_eq!(averages.len(), 3);
        assert_eq!(averages[0], 2.0); // (1+2+3)/3
        assert_eq!(averages[1], 3.0); // (2+3+4)/3
        assert_eq!(averages[2], 4.0); // (3+4+5)/3
    }

    #[test]
    fn test_correlation() {
        let x = [1.0, 2.0, 3.0, 4.0, 5.0];
        let y = [2.0, 4.0, 6.0, 8.0, 10.0];
        let corr = EnhancedSimdOps::correlation_simd(&x, &y);
        assert!((corr - 1.0).abs() < 1e-10); // Perfect positive correlation
    }

    #[test]
    fn test_log_returns() {
        let prices = [100.0, 101.0, 99.0, 102.0];
        let returns = EnhancedSimdOps::log_returns_simd(&prices);
        assert_eq!(returns.len(), 3);
        assert!((returns[0] - (101.0f64 / 100.0f64).ln()).abs() < 1e-10);
        assert!((returns[1] - (99.0f64 / 101.0f64).ln()).abs() < 1e-10);
        assert!((returns[2] - (102.0f64 / 99.0f64).ln()).abs() < 1e-10);
    }

    #[test]
    fn test_soa_processing() {
        let prices = [100.0, 101.0, 102.0, 103.0];
        let quantities = [10.0, 20.0, 30.0, 40.0];

        // Calculate notional values (price * quantity)
        let notionals =
            EnhancedSimdOps::process_order_book_levels_soa(&prices, &quantities, |price, qty| {
                price * qty
            });

        assert_eq!(notionals.len(), 4);
        assert_eq!(notionals[0], 1000.0);
        assert_eq!(notionals[1], 2020.0);
        assert_eq!(notionals[2], 3060.0);
        assert_eq!(notionals[3], 4120.0);
    }

    #[test]
    fn test_simd_aligned_buffer() {
        let mut buffer = SimdAlignedBuffer::<f64, 8>::new();
        let slice = buffer.as_mut_slice();
        slice[0] = 1.0;
        slice[7] = 8.0;

        assert_eq!(buffer.len(), 8);
        assert_eq!(buffer.as_slice()[0], 1.0);
        assert_eq!(buffer.as_slice()[7], 8.0);
    }

    #[test]
    fn test_variance_two_pass_simd_correctness() {
        // Test against known variance calculation using standard two-pass algorithm
        let values = [2.0, 4.0, 6.0, 8.0];
        let result = EnhancedSimdOps::variance_two_pass_simd(&values);

        // Calculate expected variance manually:
        // Mean = (2+4+6+8)/4 = 5.0
        // Variance = ((2-5)² + (4-5)² + (6-5)² + (8-5)²) / 4 = (9+1+1+9)/4 = 5.0
        let expected = 5.0;
        assert!(
            (result - expected).abs() < 1e-10,
            "Expected {expected}, got {result}"
        );
    }

    #[test]
    fn test_variance_two_pass_simd_aligned_vs_padded() {
        // Test with array length exactly divisible by 4 (aligned)
        let aligned_values = [1.0, 2.0, 3.0, 4.0]; // Length 4 - no padding
        let aligned_result = EnhancedSimdOps::variance_two_pass_simd(&aligned_values);

        // Test with array length requiring padding
        let padded_values = [1.0, 2.0, 3.0, 4.0, 5.0]; // Length 5 - needs padding
        let padded_result = EnhancedSimdOps::variance_two_pass_simd(&padded_values);

        // Calculate expected variances manually
        // Aligned: mean=2.5, variance=((1-2.5)²+(2-2.5)²+(3-2.5)²+(4-2.5)²)/4 = 1.25
        assert!((aligned_result - 1.25).abs() < 1e-10);

        // Padded: mean=3.0, variance=((1-3)²+(2-3)²+(3-3)²+(4-3)²+(5-3)²)/5 = 2.0
        assert!((padded_result - 2.0).abs() < 1e-10);
    }

    #[test]
    fn test_variance_two_pass_simd_edge_cases() {
        // Empty array
        let empty: &[f64] = &[];
        assert_eq!(EnhancedSimdOps::variance_two_pass_simd(empty), 0.0);

        // Single element
        let single = [42.0];
        assert_eq!(EnhancedSimdOps::variance_two_pass_simd(&single), 0.0);

        // Two identical elements
        let identical = [5.0, 5.0];
        assert_eq!(EnhancedSimdOps::variance_two_pass_simd(&identical), 0.0);

        // Two different elements
        let two_diff = [1.0, 3.0];
        // Mean = 2.0, variance = ((1-2)² + (3-2)²) / 2 = 1.0
        assert!((EnhancedSimdOps::variance_two_pass_simd(&two_diff) - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_variance_two_pass_simd_large_arrays() {
        // Test with larger arrays to verify SIMD chunking works correctly
        let large_values: Vec<f64> = (1..=100).map(|i| i as f64).collect();
        let result = EnhancedSimdOps::variance_two_pass_simd(&large_values);

        // For sequence 1..100: mean = 50.5, variance = 833.25
        let expected = 833.25;
        assert!(
            (result - expected).abs() < 1e-8,
            "Expected {expected}, got {result}"
        );

        // Test odd-length array that requires padding
        let odd_values: Vec<f64> = (1..=99).map(|i| i as f64).collect();
        let odd_result = EnhancedSimdOps::variance_two_pass_simd(&odd_values);

        // For sequence 1..99: mean = 50.0, variance = 816.6667
        let expected_odd = 816.6666666666666;
        assert!(
            (odd_result - expected_odd).abs() < 1e-8,
            "Expected {expected_odd}, got {odd_result}"
        );
    }

    #[test]
    fn test_variance_two_pass_simd_precision() {
        // Test with high precision values to ensure SIMD doesn't introduce significant errors
        let precise_values = [1.000000001, 1.000000002, 1.000000003, 1.000000004];
        let result = EnhancedSimdOps::variance_two_pass_simd(&precise_values);

        // Small variance should be calculated precisely
        assert!(result > 0.0, "Should detect small variance");
        assert!(result < 1e-15, "Small variance should be very small");

        // Test with financial data ranges (typical price values)
        let price_data = [100.25, 100.50, 99.75, 101.00, 100.125];
        let price_variance = EnhancedSimdOps::variance_two_pass_simd(&price_data);
        assert!(
            price_variance > 0.0,
            "Should calculate variance for realistic price data"
        );
    }

    #[test]
    fn test_variance_two_pass_simd_vs_welford() {
        // Compare SIMD two-pass with Welford algorithm for consistency
        let test_data = [2.5, 4.7, 1.3, 8.9, 6.1, 3.4, 7.8, 5.2];

        let simd_result = EnhancedSimdOps::variance_two_pass_simd(&test_data);
        let welford_result = EnhancedSimdOps::variance_welford_scalar(&test_data);

        // Both algorithms should give very similar results (within floating point precision)
        assert!(
            (simd_result - welford_result).abs() < 1e-12,
            "SIMD two-pass ({simd_result}) and Welford ({welford_result}) results should be nearly identical"
        );
    }
}