opendeviationbar-core 13.75.1

//! Tier 2 inter-bar features: Kyle Lambda, Burstiness, Volume moments, Garman-Klass
//!
//! Statistical features computed from lookback window.

use crate::interbar_types::TradeSnapshot;
use libm;

/// Garman-Klass volatility coefficient: 2*ln(2) - 1
/// Precomputed to avoid repeated calculation in every call
/// Exact value: 0.3862943611198906
const GARMAN_KLASS_COEFFICIENT: f64 = 0.386_294_361_119_890_6;

/// Normalize raw central moment sums into skewness and excess kurtosis.
///
/// Shared by `compute_volume_moments` and `compute_volume_moments_with_mean`
/// to eliminate duplicated normalization logic (Phase 7B dedup).
///
/// Arguments are the unnormalized sums: m2=Σ(d²), m3=Σ(d³), m4=Σ(d⁴), n=count.
#[inline]
fn normalize_moments(m2_sum: f64, m3_sum: f64, m4_sum: f64, n: f64) -> (f64, f64) {
    let n_inv = 1.0 / n;
    let m2 = m2_sum * n_inv;
    let m3 = m3_sum * n_inv;
    let m4 = m4_sum * n_inv;

    let sigma = m2.sqrt();

    if sigma < f64::EPSILON {
        return (0.0, 0.0); // All same volume
    }

    // Issue #96 Task #202: Pre-compute powers instead of powi() calls
    let sigma2 = sigma * sigma;
    let sigma3 = sigma2 * sigma;
    let sigma4 = sigma2 * sigma2;

    let skewness = m3 / sigma3;
    let kurtosis = m4 / sigma4 - 3.0; // Excess kurtosis

    (skewness, kurtosis)
}

/// Issue #96 Task #52: #[inline] for thin SIMD/scalar dispatcher (Phase 7C parity with burstiness)
#[inline]
pub fn compute_kyle_lambda(lookback: &[&TradeSnapshot]) -> f64 {
    // Issue #96 Task #148 Phase 2: Dispatch to SIMD or scalar based on feature flag
    #[cfg(feature = "simd-kyle-lambda")]
    {
        super::simd::compute_kyle_lambda_simd(lookback)
    }

    #[cfg(not(feature = "simd-kyle-lambda"))]
    {
        compute_kyle_lambda_scalar(lookback)
    }
}

/// Scalar implementation of Kyle Lambda computation (fallback/baseline).
/// Called via cfg(not(feature = "simd-kyle-lambda")) dispatcher + SIMD parity tests.
#[cfg(any(not(feature = "simd-kyle-lambda"), test, feature = "test-utils"))]
#[inline]
pub fn compute_kyle_lambda_scalar(lookback: &[&TradeSnapshot]) -> f64 {
    let n = lookback.len();

    if n < 2 {
        return 0.0;
    }

    // Issue #96 Task #210: Memoize first/last element access (scalar version)
    // Bounds guaranteed by n >= 2 check above; direct indexing avoids .first()/.last() overhead
    let first_price = lookback[0].price.to_f64();
    let last_price = lookback[n - 1].price.to_f64();

    // Adaptive computation: subsample large windows
    let (buy_vol, sell_vol) = if n > 500 {
        // For large windows (n > 500), subsample every 5th trade
        // Preserves order flow signal while reducing computation by ~5x
        lookback.iter().step_by(5).fold((0.0, 0.0), |acc, t| {
            if t.is_buyer_maker {
                (acc.0, acc.1 + t.volume.to_f64())
            } else {
                (acc.0 + t.volume.to_f64(), acc.1)
            }
        })
    } else {
        // Medium windows (10-500): Branchless ILP accumulation
        // Delegates to shared accumulate_buy_sell_branchless (Phase 4 dedup)
        super::accumulate_buy_sell_branchless(lookback)
    };

    let total_vol = buy_vol + sell_vol;
    let first_price_abs = first_price.abs();

    // Issue #96 Task #65: Coarse bounds check for extreme imbalance (early-exit optimization)
    // If one volume dominates completely (other volume ~= 0), imbalance is extreme (|imbalance| >= 1.0 - eps)
    // and we can return early without expensive normalization
    if buy_vol >= total_vol - f64::EPSILON {
        // All buys: normalized_imbalance ≈ 1.0
        return if first_price_abs > f64::EPSILON {
            (last_price - first_price) / first_price
        } else {
            0.0
        };
    } else if sell_vol >= total_vol - f64::EPSILON {
        // All sells: normalized_imbalance ≈ -1.0
        return if first_price_abs > f64::EPSILON {
            -((last_price - first_price) / first_price)
        } else {
            0.0
        };
    }

    let normalized_imbalance = if total_vol > f64::EPSILON {
        (buy_vol - sell_vol) / total_vol
    } else {
        0.0
    };

    // Issue #96 Task #208: Early-exit for zero imbalance
    // If buy_vol ≈ sell_vol (perfectly balanced), Kyle Lambda = price_change / 0 = undefined
    // Skip expensive price change calculation and return 0.0 immediately
    let imbalance_abs = normalized_imbalance.abs();
    if imbalance_abs <= f64::EPSILON {
        return 0.0; // Balanced imbalance -> Kyle Lambda = 0.0
    }

    // Issue #96 Task #203: Branchless epsilon handling using masks
    // Avoids branch misprediction penalties by checking preconditions once
    // Pattern: similar to Task #200 (OFI branchless), mask-based arithmetic
    // Branchless precondition checks: convert booleans to 0.0/1.0 masks
    let imbalance_valid = 1.0; // Already verified imbalance_abs > f64::EPSILON above
    let price_valid = if first_price_abs > f64::EPSILON {
        1.0
    } else {
        0.0
    };
    let both_valid = imbalance_valid * price_valid; // 1.0 iff both valid

    // Compute price change with guard against division by zero
    let price_change = if first_price_abs > f64::EPSILON {
        (last_price - first_price) / first_price
    } else {
        0.0
    };

    // Final result: only divide if both preconditions satisfied
    if both_valid > 0.0 {
        price_change / normalized_imbalance
    } else {
        0.0
    }
}

/// Compute Burstiness (Goh-Barabasi)
///
/// Formula: B = (sigma_tau - mu_tau) / (sigma_tau + mu_tau)
///
/// Reference: Goh & Barabasi (2008), EPL, Vol. 81, 48002
///
/// Interpretation:
/// - B = -1: Perfectly regular (periodic) arrivals
/// - B = 0: Poisson process
/// - B = +1: Maximally bursty
///
/// Issue #96 Task #52: #[inline] for thin SIMD/scalar dispatcher
#[inline]
pub fn compute_burstiness(lookback: &[&TradeSnapshot]) -> f64 {
    // Issue #96 Task #4: Dispatch to SIMD or scalar based on feature flag
    #[cfg(feature = "simd-burstiness")]
    {
        super::simd::compute_burstiness_simd(lookback)
    }

    #[cfg(not(feature = "simd-burstiness"))]
    {
        compute_burstiness_scalar(lookback)
    }
}

/// Scalar implementation of burstiness computation (fallback).
/// Uses Welford's algorithm for online variance computation (single pass, no intermediate allocation).
/// Called via cfg(not(feature = "simd-burstiness")) dispatcher + SIMD parity tests.
#[cfg(any(not(feature = "simd-burstiness"), test, feature = "test-utils"))]
#[inline]
pub fn compute_burstiness_scalar(lookback: &[&TradeSnapshot]) -> f64 {
    if lookback.len() < 2 {
        return 0.0;
    }

    // Compute mean and variance in a single pass using Welford's algorithm (Issue #96 Task #50)
    // Eliminates intermediate SmallVec allocation and second-pass variance computation
    let mut mean = 0.0;
    let mut m2 = 0.0; // Sum of squared deviations
    let mut count = 0.0;

    for i in 1..lookback.len() {
        let delta_t = (lookback[i].timestamp - lookback[i - 1].timestamp) as f64;
        count += 1.0;
        let delta = delta_t - mean;
        mean += delta / count;
        let delta2 = delta_t - mean;
        m2 += delta * delta2;
    }

    let variance = m2 / count;
    let sigma = variance.sqrt();

    // Issue #96 Task #213: Branchless epsilon check in burstiness (scalar path)
    // Eliminate branch on denominator > EPSILON by using .max() guard
    let denominator = sigma + mean;
    let numerator = sigma - mean;
    numerator / denominator.max(f64::EPSILON)
}

/// Compute volume moments (skewness and excess kurtosis)
///
/// Skewness: E[(V-mu)^3] / sigma^3 (Fisher-Pearson coefficient)
/// Excess Kurtosis: E[(V-mu)^4] / sigma^4 - 3 (normal distribution = 0)
///
/// Issue #96 Task #42: Single-pass computation avoids Vec<f64> allocation.
/// Two-phase: (1) compute mean, (2) compute moments with known mean.
#[inline]
pub fn compute_volume_moments(lookback: &[&TradeSnapshot]) -> (f64, f64) {
    let n = lookback.len() as f64;

    if n < 3.0 {
        return (0.0, 0.0);
    }

    // Issue #96: Pre-compute reciprocal — replaces 4 divisions with 1 division + 4 multiplications
    let n_inv = 1.0 / n;

    // Phase 1: Compute mean from volume stream
    let sum_vol = lookback.iter().fold(0.0, |acc, t| acc + t.volume.to_f64());
    let mu = sum_vol * n_inv;

    // Phase 2: Central moments in single pass (no Vec allocation)
    let (m2, m3, m4) = lookback.iter().fold((0.0, 0.0, 0.0), |(m2, m3, m4), t| {
        let v = t.volume.to_f64();
        let d = v - mu;
        let d2 = d * d;
        (m2 + d2, m3 + d2 * d, m4 + d2 * d2)
    });

    normalize_moments(m2, m3, m4, n)
}

/// Compute volume moments using pre-computed cache (Issue #96 Task #99)
///
/// Optimized version that reuses pre-computed volumes from LookbackCache
/// instead of converting from FixedPoint on each iteration.
/// Avoids redundant `.volume.to_f64()` calls - significant speedup when
/// computing multiple inter-bar features that need volume data.
///
/// # Performance
/// - Eliminates 2n `.volume.to_f64()` calls (Phase 1 + Phase 2 iterations)
/// - Single-pass with pre-computed data
/// - 2-5% improvement when multiple features share same lookback
#[inline]
pub fn compute_volume_moments_cached(volumes: &[f64]) -> (f64, f64) {
    let n = volumes.len() as f64;
    if n < 3.0 {
        return (0.0, 0.0);
    }
    let sum_vol: f64 = volumes.iter().sum();
    compute_volume_moments_with_mean(volumes, sum_vol / n)
}

/// Compute volume moments with pre-computed mean (Issue #96 Task #51)
///
/// Eliminates Phase 1 sum pass when total_volume is already known.
/// Called from compute_tier2_features where cache.total_volume is pre-computed.
///
/// # Performance
/// - Eliminates O(n) sum pass (0.3-0.8% speedup on 100-500 trade windows)
/// - Only central moments pass remains
#[inline]
pub fn compute_volume_moments_with_mean(volumes: &[f64], mu: f64) -> (f64, f64) {
    let n = volumes.len() as f64;

    if n < 3.0 {
        return (0.0, 0.0);
    }

    // Phase 2: Central moments in single pass
    let (m2, m3, m4) = volumes.iter().fold((0.0, 0.0, 0.0), |(m2, m3, m4), &v| {
        let d = v - mu;
        let d2 = d * d;
        (m2 + d2, m3 + d2 * d, m4 + d2 * d2)
    });

    normalize_moments(m2, m3, m4, n)
}

/// Compute Garman-Klass volatility estimator
///
/// Formula: sigma^2 = 0.5 * ln(H/L)^2 - (2*ln(2) - 1) * ln(C/O)^2
///
/// Reference: Garman & Klass (1980), Journal of Business, vol. 53, no. 1
#[inline]
pub fn compute_garman_klass(lookback: &[&TradeSnapshot]) -> f64 {
    if lookback.is_empty() {
        return 0.0;
    }

    // Issue #96 Task #210: Memoize first/last element access (Garman-Klass)
    let count = lookback.len();
    let open_price = lookback[0].price.to_f64();
    let close_price = lookback[count - 1].price.to_f64();
    let (low_price, high_price) = lookback.iter().fold((f64::MAX, f64::MIN), |acc, t| {
        let p = t.price.to_f64();
        (acc.0.min(p), acc.1.max(p))
    });

    // Delegate to OHLC variant (Phase 7A dedup)
    compute_garman_klass_with_ohlc(open_price, high_price, low_price, close_price)
}

/// Compute Garman-Klass volatility with pre-computed OHLC
///
/// Optimization: Use when OHLC data is already extracted (batch operation).
/// Avoids redundant fold operation vs compute_garman_klass().
///
/// Returns 0.0 if OHLC data is invalid.
#[inline]
pub fn compute_garman_klass_with_ohlc(open: f64, high: f64, low: f64, close: f64) -> f64 {
    // Guard: prices must be positive
    if open <= f64::EPSILON || low <= f64::EPSILON || high <= f64::EPSILON {
        return 0.0;
    }

    // Issue #96 Task #14: Use libm::log for optimized performance (1.2-1.5x speedup)
    let log_hl = libm::log(high / low);
    let log_co = libm::log(close / open);

    // Issue #96 Task #168: Optimize powi(2) to direct multiplication (0.5-1% speedup)
    let variance = 0.5 * (log_hl * log_hl) - GARMAN_KLASS_COEFFICIENT * (log_co * log_co);

    if variance > 0.0 { variance.sqrt() } else { 0.0 }
}