rusty_common/

vectorized.rs

//! Vectorization-friendly utilities for HFT calculations
//!
//! This module provides high-performance implementations of common HFT calculations
//! specifically designed to maximize compiler auto-vectorization and CPU efficiency.
//!
//! ## HFT Performance Rationale
//!
//! ### Auto-Vectorization Benefits
//! Modern CPUs can execute multiple operations simultaneously through SIMD instructions:
//! - **SSE/AVX instructions**: Process 2-8 floating-point operations per cycle
//! - **Pipeline utilization**: Keeps CPU execution units busy
//! - **Memory bandwidth**: Vectorized loads/stores maximize memory throughput
//! - **Branch reduction**: Eliminates costly conditional jumps in inner loops
//!
//! ### Critical Path Calculations
//! HFT systems require these calculations in sub-microsecond timeframes:
//! - **VWAP calculations**: Volume-weighted average price across order book levels
//! - **Volatility estimates**: Realized volatility from tick data
//! - **Order flow analysis**: Buy/sell pressure and market imbalance
//! - **Moving averages**: Exponential and simple moving averages for signals
//!
//! ## Memory Alignment Strategy
//!
//! ### Cache Line Optimization
//! ```text
//! 32-byte alignment (AVX):  ┌──── 4 x f64 ────┐
//! 64-byte alignment (Cache): ┌─────── 8 x f64 ───────┐
//! ```
//!
//! ### AlignedBuffer Types
//! - **AlignedBuffer32**: Optimized for AVX/AVX2 SIMD operations
//! - **AlignedBuffer64**: Cache-line aligned to prevent false sharing
//! - **Generic implementation**: Works with any `Copy + Default` type
//!
//! ## Compiler Vectorization Patterns
//!
//! ### Loop Structure Optimization
//! Designed to trigger auto-vectorization:
//! ```rust
//! // Vectorization-friendly pattern
//! for i in 0..values.len() {
//!     result[i] = values[i].mul_add(weights[i], accumulator);
//! }
//! ```
//!
//! ### Dependency Chain Breaking
//! Uses multiple accumulators to improve instruction-level parallelism:
//! ```rust
//! let mut sum1 = 0.0;
//! let mut sum2 = 0.0;
//! let mut sum3 = 0.0;
//! let mut sum4 = 0.0;
//!
//! // Process in chunks to break dependencies
//! for chunk in values.chunks_exact(4) {
//!     sum1 = chunk[0].mul_add(weights[0], sum1);
//!     sum2 = chunk[1].mul_add(weights[1], sum2);
//!     sum3 = chunk[2].mul_add(weights[2], sum3);
//!     sum4 = chunk[3].mul_add(weights[3], sum4);
//! }
//! ```
//!
//! ## Financial Calculation Implementations
//!
//! ### Volume-Weighted Average Price (VWAP)
//! - **Multiple accumulators**: Breaks dependency chains for better ILP
//! - **FMA operations**: Uses `mul_add` for better precision and performance
//! - **Chunk processing**: 4-element chunks optimize for SIMD width
//! - **Remainder handling**: Processes leftover elements efficiently
//!
//! ### Exponential Moving Average (EMA)
//! - **Vectorization-friendly**: Linear dependency structure
//! - **Precision optimization**: Uses `mul_add` for numerical stability
//! - **Cache efficiency**: Sequential memory access pattern
//!
//! ### Order Flow Imbalance (OFI)
//! - **Parallel processing**: Independent calculations per element
//! - **Sign operations**: Vectorizable using `signum()` function
//! - **Memory layout**: Structure-of-arrays for optimal access
//!
//! ### Realized Volatility
//! - **Squared returns**: Vectorizable multiplication operations
//! - **Multiple accumulators**: Reduces dependency chain length
//! - **Chunk processing**: Maximizes SIMD utilization
//!
//! ## Performance Characteristics
//!
//! ### Latency Improvements
//! - **2-4x speedup**: From auto-vectorization on modern CPUs
//! - **Cache efficiency**: Aligned access patterns reduce cache misses
//! - **Memory bandwidth**: Vectorized loads maximize throughput
//! - **Branch elimination**: Reduces pipeline stalls
//!
//! ### Throughput Optimization
//! - **High concurrency**: Functions are thread-safe and reentrant
//! - **Memory efficiency**: Minimal allocation through in-place operations
//! - **CPU utilization**: Keeps multiple execution units busy
//!
//! ## Integration Patterns
//!
//! ### Real-Time Trading
//! ```rust
//! // Market making: Calculate fair value
//! let fair_value = calculate_weighted_mid_price(
//!     &bid_prices, &bid_sizes,
//!     &ask_prices, &ask_sizes
//! );
//!
//! // Risk management: Estimate volatility
//! let vol = calculate_realized_volatility(&recent_returns);
//! ```
//!
//! ### Strategy Development
//! ```rust
//! // Alpha generation: Order flow analysis
//! let ofi = calculate_ofi(
//!     &bid_volumes, &ask_volumes,
//!     &bid_changes, &ask_changes
//! );
//!
//! // Signal processing: Moving averages
//! let ema = calculate_ema_vectorized(&prices, alpha);
//! ```
//!
//! ## Compiler Optimization Requirements
//!
//! ### Build Flags
//! For optimal performance, use:
//! ```toml
//! [profile.release]
//! opt-level = 3
//! lto = true
//! codegen-units = 1
//! target-cpu = "native"  # Enables AVX2/FMA if available
//! ```
//!
//! ### Target Features
//! - **AVX2**: 256-bit SIMD operations
//! - **FMA**: Fused multiply-add instructions
//! - **BMI2**: Bit manipulation for efficient indexing

/// 32-byte aligned buffer optimized for AVX/AVX2 SIMD operations
///
/// Provides guaranteed memory alignment for optimal SIMD performance in HFT calculations.
/// The 32-byte alignment matches AVX register width for efficient vectorized operations.
///
/// ## Performance Benefits
/// - **SIMD efficiency**: Aligned loads/stores avoid performance penalties
/// - **Cache optimization**: Reduces cache line splits and false sharing
/// - **Compiler hints**: Enables better auto-vectorization
/// - **Zero-cost abstraction**: Alignment overhead eliminated at compile time
///
/// ## Usage Patterns
/// ```rust
/// // Price calculations with SIMD optimization
/// let mut prices = AlignedBuffer32::<f64, 256>::new();
/// let prices_slice = prices.as_mut_slice();
/// // ... fill with market data ...
/// let vwap = calculate_vwap_vectorized(prices_slice, volumes);
/// ```
#[repr(align(32))]
pub struct AlignedBuffer32<T, const N: usize>
where
    T: Copy + Default,
{
    data: [T; N],
}

/// 64-byte aligned buffer for cache line optimization and false sharing prevention
///
/// Aligns data structures to CPU cache line boundaries (typically 64 bytes) to maximize
/// cache efficiency and prevent false sharing in multi-threaded HFT applications.
///
/// ## Cache Line Benefits
/// - **Cache efficiency**: Entire buffer fits within cache line boundaries
/// - **False sharing prevention**: Eliminates cache ping-ponging between cores
/// - **Memory bandwidth**: Maximizes utilization of cache line transfers
/// - **Prefetching optimization**: Predictable access patterns improve prefetcher efficiency
///
/// ## HFT Use Cases
/// - **Per-thread buffers**: Eliminate false sharing between trading threads
/// - **Order book snapshots**: Cache-efficient storage of market data
/// - **Statistics collections**: High-frequency counter updates
/// - **Message buffers**: Network I/O optimization
#[repr(align(64))]
pub struct AlignedBuffer64<T, const N: usize>
where
    T: Copy + Default,
{
    data: [T; N],
}

/// Legacy type alias for backward compatibility with existing f64 buffers
pub type LegacyAlignedBuffer<const N: usize> = AlignedBuffer32<f64, N>;

impl<T, const N: usize> Default for AlignedBuffer32<T, N>
where
    T: Copy + Default,
{
    fn default() -> Self {
        Self::new()
    }
}

impl<T, const N: usize> Default for AlignedBuffer64<T, N>
where
    T: Copy + Default,
{
    fn default() -> Self {
        Self::new()
    }
}

impl<T, const N: usize> AlignedBuffer32<T, N>
where
    T: Copy + Default,
{
    /// Create a new 32-byte aligned buffer with all elements initialized to default
    #[inline]
    pub fn new() -> Self {
        Self {
            data: [T::default(); N],
        }
    }

    /// Create a new 32-byte aligned buffer with all elements initialized to the given value
    #[inline]
    pub const fn with_value(value: T) -> Self {
        Self { data: [value; N] }
    }

    /// Get an immutable slice view of the buffer data
    #[inline(always)]
    pub const fn as_slice(&self) -> &[T] {
        &self.data
    }

    /// Get a mutable slice view of the buffer data
    #[inline(always)]
    pub const fn as_mut_slice(&mut self) -> &mut [T] {
        &mut self.data
    }

    /// Get the number of elements in the buffer
    #[inline(always)]
    pub const fn len(&self) -> usize {
        N
    }

    /// Check if the buffer is empty
    #[inline(always)]
    pub const fn is_empty(&self) -> bool {
        N == 0
    }

    /// Get a reference to the raw data array
    #[inline(always)]
    pub const fn data(&self) -> &[T; N] {
        &self.data
    }

    /// Get a mutable reference to the raw data array
    #[inline(always)]
    pub const fn data_mut(&mut self) -> &mut [T; N] {
        &mut self.data
    }

    /// Get the alignment of this buffer in bytes
    #[inline(always)]
    pub const fn alignment() -> usize {
        32
    }
}

impl<T, const N: usize> AlignedBuffer64<T, N>
where
    T: Copy + Default,
{
    /// Create a new 64-byte aligned buffer with all elements initialized to default
    #[inline]
    pub fn new() -> Self {
        Self {
            data: [T::default(); N],
        }
    }

    /// Create a new 64-byte aligned buffer with all elements initialized to the given value
    #[inline]
    pub const fn with_value(value: T) -> Self {
        Self { data: [value; N] }
    }

    /// Get an immutable slice view of the buffer data
    #[inline(always)]
    pub const fn as_slice(&self) -> &[T] {
        &self.data
    }

    /// Get a mutable slice view of the buffer data
    #[inline(always)]
    pub const fn as_mut_slice(&mut self) -> &mut [T] {
        &mut self.data
    }

    /// Get the number of elements in the buffer
    #[inline(always)]
    pub const fn len(&self) -> usize {
        N
    }

    /// Check if the buffer is empty
    #[inline(always)]
    pub const fn is_empty(&self) -> bool {
        N == 0
    }

    /// Get a reference to the raw data array
    #[inline(always)]
    pub const fn data(&self) -> &[T; N] {
        &self.data
    }

    /// Get a mutable reference to the raw data array
    #[inline(always)]
    pub const fn data_mut(&mut self) -> &mut [T; N] {
        &mut self.data
    }

    /// Get the alignment of this buffer in bytes
    #[inline(always)]
    pub const fn alignment() -> usize {
        64
    }
}

/// Calculate VWAP in a vectorization-friendly way
#[inline(always)]
pub fn calculate_vwap(prices: &[f64], volumes: &[f64]) -> f64 {
    debug_assert_eq!(prices.len(), volumes.len());

    // Use multiple accumulators to break dependency chains
    let mut pv_sum1 = 0.0;
    let mut pv_sum2 = 0.0;
    let mut pv_sum3 = 0.0;
    let mut pv_sum4 = 0.0;

    let mut v_sum1 = 0.0;
    let mut v_sum2 = 0.0;
    let mut v_sum3 = 0.0;
    let mut v_sum4 = 0.0;

    // Process in chunks of 4 for vectorization
    let chunks = prices.chunks_exact(4);
    let vol_chunks = volumes.chunks_exact(4);

    for (p_chunk, v_chunk) in chunks.zip(vol_chunks) {
        // Compiler will vectorize these operations
        pv_sum1 = p_chunk[0].mul_add(v_chunk[0], pv_sum1);
        pv_sum2 = p_chunk[1].mul_add(v_chunk[1], pv_sum2);
        pv_sum3 = p_chunk[2].mul_add(v_chunk[2], pv_sum3);
        pv_sum4 = p_chunk[3].mul_add(v_chunk[3], pv_sum4);

        v_sum1 += v_chunk[0];
        v_sum2 += v_chunk[1];
        v_sum3 += v_chunk[2];
        v_sum4 += v_chunk[3];
    }

    // Sum the accumulators
    let total_pv = pv_sum1 + pv_sum2 + pv_sum3 + pv_sum4;
    let total_v = v_sum1 + v_sum2 + v_sum3 + v_sum4;

    // Handle remainder
    let p_remainder = prices.chunks_exact(4).remainder();
    let v_remainder = volumes.chunks_exact(4).remainder();

    let remainder_pv: f64 = p_remainder
        .iter()
        .zip(v_remainder.iter())
        .map(|(p, v)| p * v)
        .sum();

    let remainder_v: f64 = v_remainder.iter().sum();

    (total_pv + remainder_pv) / (total_v + remainder_v)
}

/// Calculate exponential moving average with vectorization
#[inline(always)]
pub fn calculate_ema_vectorized(data: &[f64], alpha: f64) -> Vec<f64> {
    if data.is_empty() {
        return Vec::new();
    }

    let mut result = Vec::with_capacity(data.len());
    result.push(data[0]);

    let one_minus_alpha = 1.0 - alpha;

    // This loop pattern allows for better vectorization
    for i in 1..data.len() {
        let ema = alpha.mul_add(data[i], one_minus_alpha * result[i - 1]);
        result.push(ema);
    }

    result
}

/// Calculate order flow imbalance (OFI) with vectorization
#[inline(always)]
pub fn calculate_ofi(
    bid_volumes: &[f64],
    ask_volumes: &[f64],
    bid_price_changes: &[f64],
    ask_price_changes: &[f64],
) -> Vec<f64> {
    debug_assert_eq!(bid_volumes.len(), ask_volumes.len());
    debug_assert_eq!(bid_volumes.len(), bid_price_changes.len());
    debug_assert_eq!(bid_volumes.len(), ask_price_changes.len());

    let mut ofi = vec![0.0; bid_volumes.len()];

    // Vectorizable main loop
    for i in 0..bid_volumes.len() {
        let bid_component = bid_volumes[i] * bid_price_changes[i].signum();
        let ask_component = ask_volumes[i] * ask_price_changes[i].signum();
        ofi[i] = bid_component - ask_component;
    }

    ofi
}

/// Calculate realized volatility with vectorization
#[inline(always)]
pub fn calculate_realized_volatility(returns: &[f64]) -> f64 {
    if returns.is_empty() {
        return 0.0;
    }

    // Use multiple accumulators for better pipelining
    let mut sum1 = 0.0;
    let mut sum2 = 0.0;
    let mut sum3 = 0.0;
    let mut sum4 = 0.0;

    let chunks = returns.chunks_exact(4);
    let remainder = chunks.remainder();

    for chunk in chunks {
        // Squared returns - vectorizable
        sum1 = chunk[0].mul_add(chunk[0], sum1);
        sum2 = chunk[1].mul_add(chunk[1], sum2);
        sum3 = chunk[2].mul_add(chunk[2], sum3);
        sum4 = chunk[3].mul_add(chunk[3], sum4);
    }

    // Handle remainder
    let remainder_sum: f64 = remainder.iter().map(|r| r * r).sum();

    let total_sum = sum1 + sum2 + sum3 + sum4 + remainder_sum;
    (total_sum / returns.len() as f64).sqrt()
}

/// Calculate weighted mid price with vectorization
#[inline(always)]
pub fn calculate_weighted_mid_price(
    bid_prices: &[f64],
    bid_sizes: &[f64],
    ask_prices: &[f64],
    ask_sizes: &[f64],
) -> f64 {
    debug_assert_eq!(bid_prices.len(), bid_sizes.len());
    debug_assert_eq!(ask_prices.len(), ask_sizes.len());

    // Calculate total bid volume
    let total_bid_volume: f64 = bid_sizes.iter().sum();
    let total_ask_volume: f64 = ask_sizes.iter().sum();

    if total_bid_volume == 0.0 || total_ask_volume == 0.0 {
        return 0.0;
    }

    // Weighted bid price
    let weighted_bid = calculate_weighted_average(bid_prices, bid_sizes);

    // Weighted ask price
    let weighted_ask = calculate_weighted_average(ask_prices, ask_sizes);

    // Volume-weighted mid price
    let bid_weight = total_bid_volume / (total_bid_volume + total_ask_volume);
    let ask_weight = total_ask_volume / (total_bid_volume + total_ask_volume);

    weighted_bid * bid_weight + weighted_ask * ask_weight
}

/// Helper function for weighted average calculation
#[inline(always)]
fn calculate_weighted_average(values: &[f64], weights: &[f64]) -> f64 {
    debug_assert_eq!(values.len(), weights.len());

    let mut weighted_sum = 0.0;
    let mut weight_sum = 0.0;

    // Vectorizable loop
    for i in 0..values.len() {
        weighted_sum = values[i].mul_add(weights[i], weighted_sum);
        weight_sum += weights[i];
    }

    weighted_sum / weight_sum
}

/// Calculate order book imbalance with vectorization
#[inline(always)]
pub fn calculate_order_book_imbalance(bid_volumes: &[f64], ask_volumes: &[f64]) -> f64 {
    let total_bid: f64 = bid_volumes.iter().sum();
    let total_ask: f64 = ask_volumes.iter().sum();

    if total_bid + total_ask == 0.0 {
        return 0.0;
    }

    (total_bid - total_ask) / (total_bid + total_ask)
}

/// Fast moving sum calculation
#[inline(always)]
pub fn moving_sum(data: &[f64], window: usize) -> Vec<f64> {
    if data.len() < window || window == 0 {
        return vec![];
    }

    let mut result = Vec::with_capacity(data.len() - window + 1);

    // Initial sum - vectorizable
    let mut sum: f64 = data[..window].iter().sum();
    result.push(sum);

    // Rolling window - efficient single pass
    for i in window..data.len() {
        sum += data[i] - data[i - window];
        result.push(sum);
    }

    result
}

/// Vectorized returns calculation
#[inline(always)]
pub fn calculate_returns(prices: &[f64]) -> Vec<f64> {
    if prices.len() < 2 {
        return vec![];
    }

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

    // Vectorizable loop for returns
    for i in 1..prices.len() {
        let return_val = (prices[i] - prices[i - 1]) / prices[i - 1];
        returns.push(return_val);
    }

    returns
}

/// Vectorized log returns calculation
#[inline(always)]
pub fn calculate_log_returns(prices: &[f64]) -> Vec<f64> {
    if prices.len() < 2 {
        return vec![];
    }

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

    // Vectorizable loop for log returns
    for i in 1..prices.len() {
        let log_return = (prices[i] / prices[i - 1]).ln();
        returns.push(log_return);
    }

    returns
}

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

    #[test]
    fn test_calculate_vwap() {
        let prices = vec![100.0, 101.0, 102.0, 103.0, 104.0];
        let volumes = vec![1000.0, 1500.0, 2000.0, 1200.0, 1800.0];

        let vwap = calculate_vwap(&prices, &volumes);

        // Manual calculation for verification
        let total_pv: f64 = prices.iter().zip(volumes.iter()).map(|(p, v)| p * v).sum();
        let total_v: f64 = volumes.iter().sum();
        let expected = total_pv / total_v;

        eprintln!("VWAP calculated: {vwap}");
        eprintln!("VWAP expected: {expected}");
        eprintln!("Difference: {}", (vwap - expected).abs());

        // The manual calculation gives us ~102.1733
        assert!((vwap - 102.173333).abs() < 0.01);
    }

    #[test]
    fn test_calculate_ofi() {
        let bid_vols = vec![100.0, 150.0, 200.0];
        let ask_vols = vec![120.0, 130.0, 180.0];
        let bid_changes = vec![0.1, -0.1, 0.2];
        let ask_changes = vec![-0.1, 0.1, -0.2];

        let ofi = calculate_ofi(&bid_vols, &ask_vols, &bid_changes, &ask_changes);
        assert_eq!(ofi.len(), 3);
    }

    #[test]
    fn test_realized_volatility() {
        let returns = vec![0.01, -0.02, 0.015, -0.01, 0.02];
        let vol = calculate_realized_volatility(&returns);
        assert!(vol > 0.0);
    }
}