//! 42-Dimension Feature Extraction for ML Models
//!
//! This module implements streamlined feature engineering for HFT ML models,
//! extracting 42 features per OHLCV bar:
//! - 42 market features (40 base + 2 regime)
//! - 2 regime features (ADX trend strength, CUSUM direction) for GPU-native regime classification
//!
//! ## Performance
//! - Target: <1ms per bar for 42 features
//! - Memory: ~336 bytes per bar (42 × f64)
//! - Uses rolling windows (VecDeque) for O(1) amortized complexity
//!
//! ## Architecture
//! ```rust
//! use ml::features::extraction::extract_ml_features;
//! use ml::data_loader::RealDataLoader;
//!
//! let loader = RealDataLoader::new();
//! let bars = loader.load_ohlcv_bars("ES.FUT").await?;
//! let features = extract_ml_features(&bars)?;  // Vec<[f64; 42]>
//! ```

use crate::features::microstructure::{AmihudIlliquidity, CorwinSchultzSpread, RollMeasure};
pub use crate::types::OHLCVBar;
use anyhow::{Context, Result};
use chrono::{Datelike, Timelike};
use common::features::{BollingerBands, ATR, EMA, MACD, RSI};

// WAVE 8 AGENT 37: Import Wave D feature modules
use crate::features::regime_adaptive::RegimeAdaptiveFeatures;
use crate::features::regime_adx::RegimeADXFeatures;
use crate::features::regime_cusum::RegimeCUSUMFeatures;
use crate::features::regime_transition::RegimeTransitionFeatures;
use std::collections::VecDeque;

// WAVE 4: Import OFI calculator for TRUE OFI features (8 features from MBP-10 data)
use crate::features::ofi_calculator::OFICalculator;
use data::providers::databento::mbp10::Mbp10Snapshot;

/// Feature extraction version for cache invalidation
///
/// Increment this version string ONLY when the feature extraction logic changes
/// in a way that affects the output. This is used by the caching system to
/// automatically invalidate stale caches.
///
/// Version History:
/// - 1.0.0: Initial 51-feature implementation (43 base + 8 OFI)
/// - 2.0.0: Remove dead features (3 regime + 8 OFI) → 42 market features (40 base + 2 regime: ADX, CUSUM)
/// - 3.0.0: Add ADX + CUSUM direction for GPU-native regime classification → 42 market features
pub const FEATURE_EXTRACTION_VERSION: &str = "3.0.0";

/// Feature extraction result: 42-dimensional feature vector per bar
///
/// Layout: OHLCV(5) + Technical(5) + Patterns(6) + Volume(6) + Time(5) + Statistical(13) + Regime(2) = 42
pub type FeatureVector = [f64; 42];

/// Main feature extraction function: Converts OHLCV bars to 42-dim feature vectors
///
/// ## Arguments
/// - `bars`: Input OHLCV bars from real data loader
///
/// ## Returns
/// - `Vec<FeatureVector>`: 42-dim feature vectors per bar (after warmup period)
///
/// ## Feature Breakdown (42 features)
/// - Features 0-4: OHLCV (5 features)
/// - Features 5-9: Technical indicators (5 features)
/// - Features 10-15: Price patterns (6 features)
/// - Features 16-21: Volume features (6 features)
/// - Features 22-26: Time-based features (5 features)
/// - Features 27-39: Statistical features (13 features)
/// - Feature 40: ADX(14) normalized [0,1] (regime trend strength)
/// - Feature 41: CUSUM direction [-1,1] (regime change direction)
///
/// ## Warmup Period
/// Requires minimum 50 bars for rolling windows. Returns feature vectors only after warmup.
pub fn extract_ml_features(bars: &[OHLCVBar]) -> Result<Vec<FeatureVector>> {
    if bars.is_empty() {
        anyhow::bail!("Cannot extract features from empty bar sequence");
    }

    const WARMUP_PERIOD: usize = 50;
    if bars.len() < WARMUP_PERIOD {
        anyhow::bail!(
            "Insufficient data: {} bars provided, {} required for warmup",
            bars.len(),
            WARMUP_PERIOD
        );
    }

    let mut extractor = FeatureExtractor::new();
    let mut feature_vectors = Vec::with_capacity(bars.len() - WARMUP_PERIOD);

    // Feed bars sequentially to build rolling windows
    for (i, bar) in bars.iter().enumerate() {
        extractor.update(bar)?;

        // Start extracting features after warmup
        if i >= WARMUP_PERIOD {
            feature_vectors.push(extractor.extract_current_features()?);
        }
    }

    Ok(feature_vectors)
}

/// Stateful feature extractor with rolling windows for O(1) amortized complexity
// WAVE 7 AGENT 29C: Made public to allow custom Wave C extraction in DQN trainer
#[derive(Debug)]
pub struct FeatureExtractor {
    /// Rolling window of bars (max 260 for 52-week approximation)
    bars: VecDeque<OHLCVBar>,
    /// Technical indicator calculator (reuse from ml_training_service)
    indicators: TechnicalIndicatorState,
    /// Roll Measure (effective spread estimator)
    roll_measure: RollMeasure,
    /// Amihud Illiquidity (price impact measure)
    amihud_illiquidity: AmihudIlliquidity,
    /// Corwin-Schultz Spread (high-low volatility decomposition)
    corwin_schultz_spread: CorwinSchultzSpread,

    // WAVE 8 AGENT 37: Wave D feature extractors (indices 201-224, 24 features)
    /// CUSUM regime detection features (indices 201-210, 10 features)
    regime_cusum: RegimeCUSUMFeatures,
    /// ADX directional indicators (indices 211-215, 5 features)
    regime_adx: RegimeADXFeatures,
    /// Transition probabilities (indices 216-220, 5 features)
    regime_transition: RegimeTransitionFeatures,
    /// Adaptive position/stop-loss metrics (indices 221-224, 4 features)
    regime_adaptive: RegimeAdaptiveFeatures,

    // WAVE 4: OFI calculator for TRUE OFI features (indices 46-53, 8 features)
    /// Order Flow Imbalance calculator (stateful, requires MBP-10 data)
    ofi_calculator: OFICalculator,

    /// Last ADX(14) value from regime_adx extractor, range [0, 100]
    last_adx: f64,
    /// Last CUSUM direction from regime_cusum extractor, range [-1, 1]
    last_cusum_direction: f64,
}

impl FeatureExtractor {
    pub fn new() -> Self {
        Self {
            bars: VecDeque::with_capacity(260),
            indicators: TechnicalIndicatorState::new(),
            roll_measure: RollMeasure::new(),
            amihud_illiquidity: AmihudIlliquidity::default(),
            corwin_schultz_spread: CorwinSchultzSpread::new(),
            // WAVE 8 AGENT 37: Initialize Wave D extractors
            regime_cusum: RegimeCUSUMFeatures::new(0.0, 1.0, 0.5, 4.0),
            regime_adx: RegimeADXFeatures::new(14),
            regime_transition: RegimeTransitionFeatures::new(4, 0.1),
            regime_adaptive: RegimeAdaptiveFeatures::new(20, 100_000.0, 14),
            // WAVE 4: Initialize OFI calculator
            ofi_calculator: OFICalculator::new(),
            // Regime feature cache (populated by update(), consumed by extract_current_features_v2())
            last_adx: 0.0,
            last_cusum_direction: 0.0,
        }
    }

    pub fn update(&mut self, bar: &OHLCVBar) -> Result<()> {
        // Add to rolling window
        self.bars.push_back(bar.clone());
        if self.bars.len() > 260 {
            self.bars.pop_front();
        }

        // Update technical indicators
        self.indicators.update(bar)?;

        // Update microstructure features
        self.roll_measure.update(bar.close);
        self.amihud_illiquidity.update(bar.close, bar.volume);
        self.corwin_schultz_spread
            .update(bar.high, bar.low, bar.close);

        // Update regime feature extractors (ADX + CUSUM for GPU-native regime classification)
        let adx_features = self.regime_adx.update(bar);
        self.last_adx = adx_features[0]; // ADX(14) value, range [0, 100]

        let close_return = if self.bars.len() > 1 {
            let prev_close = self.bars[self.bars.len() - 2].close;
            safe_log_return(bar.close, prev_close)
        } else {
            0.0
        };
        let cusum_features = self.regime_cusum.update(close_return);
        self.last_cusum_direction = cusum_features[3]; // Direction: -1, 0, or +1

        Ok(())
    }

    /// Extract 42 market features for the current bar state.
    ///
    /// Returns FeatureVector ([f64; 42]) containing:
    /// - OHLCV (5), Technical (5), Price Patterns (6), Volume (6)
    /// - Time (5), Statistical (13), Regime (2: ADX, CUSUM direction)
    pub fn extract_current_features_v2(&mut self) -> Result<FeatureVector> {
        let mut features = [0.0; 42];
        let mut offset = 0;

        // OHLCV (5 features, indices 0-4)
        self.extract_ohlcv_features(&mut features[offset..offset+5])?;
        offset += 5;

        // Technical (5 features, indices 5-9)
        self.extract_technical_features_v2(&mut features[offset..offset+5])?;
        offset += 5;

        // Price Patterns (6 features, indices 10-15)
        self.extract_price_patterns_v2(&mut features[offset..offset+6])?;
        offset += 6;

        // Volume (6 features, indices 16-21)
        self.extract_volume_features_v2(&mut features[offset..offset+6])?;
        offset += 6;

        // Time (5 features, indices 22-26)
        self.extract_time_features_v2(&mut features[offset..offset+5])?;
        offset += 5;

        // Statistical (13 features, indices 27-39)
        self.extract_statistical_features_v2(&mut features[offset..offset+13])?;
        offset += 13;

        // Regime (2 features, indices 40-41)
        features[offset] = safe_normalize(self.last_adx, 0.0, 100.0); // ADX [0,1]
        offset += 1;
        features[offset] = self.last_cusum_direction; // CUSUM direction [-1, 1]
        offset += 1;

        debug_assert_eq!(offset, 42, "Expected 42 core features, got {}", offset);

        // Validate no NaN/Inf
        self.validate_features(&features)?;

        Ok(features)
    }

    /// Extract 42 market features (OFI features removed in Task 5.1).
    ///
    /// MBP-10 snapshots are accepted for API compatibility but ignored.
    /// OFI features were always zero-padded when MBP-10 data was unavailable
    /// and have been removed as dead features.
    ///
    /// # Arguments
    /// * `mbp10_snapshots` - Ignored (kept for API compatibility)
    ///
    /// # Returns
    /// * `Result<FeatureVector>` - 42-element feature vector
    pub fn extract_current_features_with_ofi(
        &mut self,
        _mbp10_snapshots: &[Mbp10Snapshot],
    ) -> Result<FeatureVector> {
        // OFI features removed — just return the 42 market features
        self.extract_current_features_v2()
    }


    /// Extract technical indicators (5) - Version 2 (REDUCED from 10)
    ///
    /// Features 5-9:
    /// - Index 5: RSI(14) normalized [0,1]
    /// - Index 6: MACD histogram only (captures divergence)
    /// - Index 7: Bollinger upper band
    /// - Index 8: Bollinger lower band
    /// - Index 9: ATR(14)
    fn extract_technical_features_v2(&self, out: &mut [f64]) -> Result<()> {
        let indicators = &self.indicators;

        out[0] = safe_normalize(indicators.last_rsi, 0.0, 100.0); // RSI
        out[1] = safe_clip(indicators.last_macd.2, -3.0, 3.0); // MACD histogram only
        out[2] = safe_clip(indicators.last_bollinger.1, -3.0, 3.0); // BB upper
        out[3] = safe_clip(indicators.last_bollinger.2, -3.0, 3.0); // BB lower
        out[4] = safe_normalize(indicators.last_atr.max(0.001).ln(), -7.0, 9.0); // ATR log-normalized (ln(0.001)≈-6.9, ln(5000)≈8.5)

        Ok(())
    }

    /// Extract price patterns (6) - Version 2 (REDUCED from 60)
    ///
    /// Features 10-15:
    /// - Index 10: Simple return (close to prev close)
    /// - Index 11: Intraday return (close to open)
    /// - Index 12: Overnight return (open to prev close)
    /// - Index 13: Close to SMA(20) ratio
    /// - Index 14: Close to SMA(50) ratio
    /// - Index 15: Linear regression slope(20)
    fn extract_price_patterns_v2(&self, out: &mut [f64]) -> Result<()> {
        let bar = self.bars.back().context("No current bar")?;
        let mut idx = 0;

        // Returns (3)
        if self.bars.len() > 1 {
            let prev = &self.bars[self.bars.len() - 2];
            out[idx] = safe_log_return(bar.close, prev.close); // Simple return
            idx += 1;
            out[idx] = safe_log_return(bar.close, bar.open); // Intraday return
            idx += 1;
            out[idx] = safe_log_return(bar.open, prev.close); // Overnight return
            idx += 1;
        } else {
            idx += 3;
        }

        // Moving average ratios (2)
        out[idx] = if self.bars.len() >= 20 {
            let sma = self.compute_sma(20);
            safe_clip((bar.close / sma) - 1.0, -0.5, 0.5)
        } else {
            0.0
        };
        idx += 1;

        out[idx] = if self.bars.len() >= 50 {
            let sma = self.compute_sma(50);
            safe_clip((bar.close / sma) - 1.0, -0.5, 0.5)
        } else {
            0.0
        };
        idx += 1;

        // Linear regression slope (1)
        out[idx] = if self.bars.len() >= 20 {
            let slope = self.compute_linear_regression_slope(20);
            safe_clip(slope, -0.1, 0.1)
        } else {
            0.0
        };

        Ok(())
    }

    /// Extract volume features (6) - Version 2 (REDUCED from 40)
    ///
    /// Features 16-21:
    /// - Index 16: Volume ratio to SMA(20)
    /// - Index 17: Volume spike indicator (>2x avg)
    /// - Index 18: VWAP(20)
    /// - Index 19: VWAP deviation
    /// - Index 20: Price-volume product
    /// - Index 21: Price-volume correlation(20)
    fn extract_volume_features_v2(&self, out: &mut [f64]) -> Result<()> {
        let bar = self.bars.back().context("No current bar")?;
        let mut idx = 0;

        // Volume ratio to SMA(20)
        out[idx] = if self.bars.len() >= 20 {
            let vol_sma = self.compute_volume_sma(20);
            safe_clip((bar.volume / vol_sma) - 1.0, -2.0, 2.0)
        } else {
            0.0
        };
        idx += 1;

        // Volume spike indicator
        out[idx] = if self.bars.len() >= 5 {
            let avg_vol = self.compute_volume_sma(5);
            if bar.volume > avg_vol * 2.0 {
                1.0
            } else {
                0.0
            }
        } else {
            0.0
        };
        idx += 1;

        // VWAP(20) - normalized as ratio to current close
        out[idx] = if self.bars.len() >= 20 && bar.close > 0.0 {
            let vwap = self.compute_vwap(20);
            safe_clip((vwap / bar.close) - 1.0, -0.1, 0.1)
        } else {
            0.0
        };
        idx += 1;

        // VWAP deviation
        out[idx] = safe_clip(
            (bar.close / (self.compute_vwap(20) + 1e-8)) - 1.0,
            -0.1,
            0.1,
        );
        idx += 1;

        // Price-volume product
        out[idx] = if self.bars.len() > 1 {
            let ret = safe_log_return(bar.close, self.bars[self.bars.len() - 2].close);
            ret * safe_normalize(bar.volume.max(1.0).ln(), 0.0, 18.0)
        } else {
            0.0
        };
        idx += 1;

        // Price-volume correlation(20)
        out[idx] = self.compute_price_volume_correlation(20);

        Ok(())
    }

    /// Extract time-based features (5) - Version 2 (REDUCED from 10)
    ///
    /// Features 22-26:
    /// - Index 22: Hour of day (0-23 normalized)
    /// - Index 23: Day of week (0-6 normalized)
    /// - Index 24: Is market open (binary)
    /// - Index 25: Minutes since open
    /// - Index 26: Minutes to close
    fn extract_time_features_v2(&self, out: &mut [f64]) -> Result<()> {
        let bar = self.bars.back().context("No current bar")?;
        let dt = bar.timestamp;

        out[0] = safe_normalize(dt.hour() as f64, 0.0, 23.0);
        out[1] = safe_normalize(dt.weekday().num_days_from_monday() as f64, 0.0, 6.0);
        out[2] = if dt.hour() >= 9 && dt.hour() < 16 {
            1.0
        } else {
            0.0
        };
        out[3] = safe_normalize(
            (dt.hour() as f64 - 9.0) * 60.0 + dt.minute() as f64,
            0.0,
            420.0,
        );
        out[4] = safe_normalize(
            (16.0 - dt.hour() as f64) * 60.0 - dt.minute() as f64,
            0.0,
            420.0,
        );

        Ok(())
    }

    /// Extract statistical features (13) - Version 2 (REDUCED from 26)
    ///
    /// Features 27-39:
    /// - Indices 27-28: Z-score(10), Z-score(20)
    /// - Indices 29-30: Percentile(10), Percentile(20)
    /// - Indices 31-33: Autocorr lag-1, lag-5, lag-10
    /// - Indices 34-36: Skewness(5), Skewness(10), Skewness(20)
    /// - Indices 37-39: Kurtosis(5), Kurtosis(10), Kurtosis(20)
    fn extract_statistical_features_v2(&self, out: &mut [f64]) -> Result<()> {
        let bar = self.bars.back().context("No current bar")?;
        let mut idx = 0;

        // Z-scores (2)
        for period in [10, 20] {
            if self.bars.len() >= period {
                let mean = self.compute_sma(period);
                let std = self.compute_std(period);
                out[idx] = safe_clip((bar.close - mean) / (std + 1e-8), -3.0, 3.0);
                idx += 1;
            } else {
                idx += 1;
            }
        }

        // Percentile ranks (2)
        for period in [10, 20] {
            if self.bars.len() >= period {
                let min = self.compute_min(period);
                let max = self.compute_max(period);
                out[idx] = safe_clip((bar.close - min) / (max - min + 1e-8), 0.0, 1.0);
                idx += 1;
            } else {
                idx += 1;
            }
        }

        // Autocorrelations (3)
        for lag in [1, 5, 10] {
            out[idx] = if self.bars.len() > lag {
                self.compute_autocorr(lag)
            } else {
                0.0
            };
            idx += 1;
        }

        // Skewness (3)
        for period in [5, 10, 20] {
            out[idx] = self.compute_skewness(period);
            idx += 1;
        }

        // Kurtosis (3)
        for period in [5, 10, 20] {
            out[idx] = self.compute_kurtosis(period);
            idx += 1;
        }

        Ok(())
    }

    // extract_regime_features_v2 removed (Task 5.1): always returned zeros

    /// Extract current 42-dim market features.
    ///
    /// Returns the same result as extract_current_features_v2().
    pub fn extract_current_features(&mut self) -> Result<FeatureVector> {
        self.extract_current_features_v2()
    }

    /// Extract OHLCV features (5): Normalized raw price/volume data
    fn extract_ohlcv_features(&self, out: &mut [f64]) -> Result<()> {
        let bar = self.bars.back().context("No current bar")?;

        // Normalize using log returns and volume ratio
        let prev_close = if self.bars.len() > 1 {
            self.bars[self.bars.len() - 2].close
        } else {
            bar.close
        };

        out[0] = safe_log_return(bar.open, prev_close); // Open relative to prev close
        out[1] = safe_log_return(bar.high, prev_close); // High relative to prev close
        out[2] = safe_log_return(bar.low, prev_close); // Low relative to prev close
        out[3] = safe_log_return(bar.close, prev_close); // Close return
        out[4] = safe_normalize(bar.volume.max(1.0).ln(), 0.0, 18.0); // Volume log-normalized (ln(1)=0, ln(10M)≈16.1)

        Ok(())
    }

    /// Validate no NaN/Inf in feature vector
    fn validate_features(&self, features: &[f64]) -> Result<()> {
        for (i, &val) in features.iter().enumerate() {
            if !val.is_finite() {
                anyhow::bail!("Invalid feature at index {}: {}", i, val);
            }
        }
        Ok(())
    }

    // ===== Helper Methods =====

    fn compute_sma(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        let sum: f64 = self.bars.iter().skip(start).map(|b| b.close).sum();
        sum / period as f64
    }

    fn compute_std(&self, period: usize) -> f64 {
        let mean = self.compute_sma(period);
        let start = self.bars.len().saturating_sub(period);
        let variance: f64 = self
            .bars
            .iter()
            .skip(start)
            .map(|b| (b.close - mean).powi(2))
            .sum::<f64>()
            / period as f64;
        variance.sqrt()
    }

    fn compute_min(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        self.bars
            .iter()
            .skip(start)
            .map(|b| b.close)
            .fold(f64::INFINITY, f64::min)
    }

    fn compute_max(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        self.bars
            .iter()
            .skip(start)
            .map(|b| b.close)
            .fold(f64::NEG_INFINITY, f64::max)
    }

    fn compute_median(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        let mut values: Vec<f64> = self.bars.iter().skip(start).map(|b| b.close).collect();
        values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        values[values.len() / 2]
    }

    fn compute_volume_sma(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        let sum: f64 = self.bars.iter().skip(start).map(|b| b.volume).sum();
        sum / period as f64
    }

    fn compute_volume_std(&self, period: usize) -> f64 {
        let mean = self.compute_volume_sma(period);
        let start = self.bars.len().saturating_sub(period);
        let variance: f64 = self
            .bars
            .iter()
            .skip(start)
            .map(|b| (b.volume - mean).powi(2))
            .sum::<f64>()
            / period as f64;
        variance.sqrt()
    }

    fn compute_vwap(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        let (weighted_sum, volume_sum): (f64, f64) = self
            .bars
            .iter()
            .skip(start)
            .map(|b| (b.close * b.volume, b.volume))
            .fold((0.0, 0.0), |(ws, vs), (w, v)| (ws + w, vs + v));
        weighted_sum / (volume_sum + 1e-8)
    }

    fn compute_momentum(&self, period: usize) -> f64 {
        if self.bars.len() > period {
            let curr = match self.bars.back() {
                Some(b) => b.close,
                None => return 0.0,
            };
            let prev = self.bars[self.bars.len() - period - 1].close;
            if prev <= 0.0 || !prev.is_finite() {
                return 0.0;
            }
            let momentum = (curr - prev) / prev;
            if !momentum.is_finite() {
                return 0.0;
            }
            momentum
        } else {
            0.0
        }
    }

    fn compute_linear_regression_slope(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let start = self.bars.len() - period;
        let n = period as f64;
        let sum_x = (n * (n - 1.0)) / 2.0; // 0 + 1 + ... + (n-1)
        let sum_x2 = (n * (n - 1.0) * (2.0 * n - 1.0)) / 6.0; // Sum of squares
        let mut sum_y = 0.0;
        let mut sum_xy = 0.0;
        for (i, bar) in self.bars.iter().skip(start).enumerate() {
            sum_y += bar.close;
            sum_xy += i as f64 * bar.close;
        }
        (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x * sum_x)
    }

    fn compute_autocorr(&self, lag: usize) -> f64 {
        if self.bars.len() <= lag {
            return 0.0;
        }
        let n = self.bars.len() - lag;
        let mean: f64 = self.bars.iter().map(|b| b.close).sum::<f64>() / self.bars.len() as f64;
        let mut numerator = 0.0;
        let mut denominator = 0.0;
        for i in 0..n {
            numerator += (self.bars[i].close - mean) * (self.bars[i + lag].close - mean);
        }
        for bar in self.bars.iter() {
            denominator += (bar.close - mean).powi(2);
        }
        numerator / (denominator + 1e-8)
    }

    // Missing helper methods implementation
    fn compute_distance_to_high(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let max = self.compute_max(period);
        let current = match self.bars.back() {
            Some(b) => b.close,
            None => return 0.0,
        };
        safe_clip((current - max) / current, -0.5, 0.0)
    }

    fn compute_distance_to_low(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let min = self.compute_min(period);
        let current = match self.bars.back() {
            Some(b) => b.close,
            None => return 0.0,
        };
        safe_clip((current - min) / current, 0.0, 0.5)
    }

    fn compute_percentile_rank(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.5;
        }
        let current = match self.bars.back() {
            Some(b) => b.close,
            None => return 0.5,
        };
        let start = self.bars.len().saturating_sub(period);
        let count_below = self
            .bars
            .iter()
            .skip(start)
            .filter(|b| b.close < current)
            .count();
        count_below as f64 / period as f64
    }

    fn compute_consecutive_highs(&self) -> f64 {
        let mut count = 0;
        if self.bars.len() < 2 {
            return 0.0;
        }
        for i in (0..self.bars.len() - 1).rev() {
            if self.bars[i + 1].close > self.bars[i].close {
                count += 1;
            } else {
                break;
            }
        }
        safe_normalize(count as f64, 0.0, 10.0)
    }

    fn compute_consecutive_lows(&self) -> f64 {
        let mut count = 0;
        if self.bars.len() < 2 {
            return 0.0;
        }
        for i in (0..self.bars.len() - 1).rev() {
            if self.bars[i + 1].close < self.bars[i].close {
                count += 1;
            } else {
                break;
            }
        }
        safe_normalize(count as f64, 0.0, 10.0)
    }

    fn compute_trend_quality(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let slope = self.compute_linear_regression_slope(period);
        let std = self.compute_std(period);
        let mean = self.compute_sma(period);
        safe_clip(slope.abs() / (std / mean + 1e-8), 0.0, 1.0)
    }

    fn compute_roc(&self, period: usize) -> f64 {
        if self.bars.len() <= period {
            return 0.0;
        }
        let current = match self.bars.back() {
            Some(b) => b.close,
            None => return 0.0,
        };
        let prev = self.bars[self.bars.len() - period - 1].close;
        if prev <= 0.0 || !prev.is_finite() {
            return 0.0;
        }
        let roc = (current - prev) / prev;
        safe_clip(roc, -0.5, 0.5)
    }

    fn compute_price_acceleration(&self) -> f64 {
        if self.bars.len() < 3 {
            return 0.0;
        }
        let curr = match self.bars.back() {
            Some(b) => b.close,
            None => return 0.0,
        };
        let prev1 = self.bars[self.bars.len() - 2].close;
        let prev2 = self.bars[self.bars.len() - 3].close;
        let vel1 = curr - prev1;
        let vel2 = prev1 - prev2;
        safe_clip(vel1 - vel2, -1.0, 1.0)
    }

    fn compute_price_velocity(&self) -> f64 {
        if self.bars.len() < 2 {
            return 0.0;
        }
        let curr = match self.bars.back() {
            Some(b) => b.close,
            None => return 0.0,
        };
        let prev = self.bars[self.bars.len() - 2].close;
        safe_clip(curr - prev, -1.0, 1.0)
    }

    fn compute_body_ratio(&self) -> f64 {
        let bar = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let body = (bar.close - bar.open).abs();
        let range = bar.high - bar.low + 1e-8;
        safe_clip(body / range, 0.0, 1.0)
    }

    fn compute_upper_shadow_ratio(&self) -> f64 {
        let bar = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let upper_shadow = bar.high - bar.close.max(bar.open);
        let range = bar.high - bar.low + 1e-8;
        safe_clip(upper_shadow / range, 0.0, 1.0)
    }

    fn compute_lower_shadow_ratio(&self) -> f64 {
        let bar = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let lower_shadow = bar.close.min(bar.open) - bar.low;
        let range = bar.high - bar.low + 1e-8;
        safe_clip(lower_shadow / range, 0.0, 1.0)
    }

    fn compute_doji_indicator(&self) -> f64 {
        let bar = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let body = (bar.close - bar.open).abs();
        let range = bar.high - bar.low + 1e-8;
        if body / range < 0.1 {
            1.0
        } else {
            0.0
        }
    }

    fn compute_hammer_indicator(&self) -> f64 {
        let bar = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let body = (bar.close - bar.open).abs();
        let lower_shadow = bar.close.min(bar.open) - bar.low;
        let range = bar.high - bar.low + 1e-8;
        if lower_shadow > body * 2.0 && body / range > 0.1 {
            1.0
        } else {
            0.0
        }
    }

    fn compute_engulfing_indicator(&self) -> f64 {
        if self.bars.len() < 2 {
            return 0.0;
        }
        let curr = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let prev = &self.bars[self.bars.len() - 2];
        let curr_body = (curr.close - curr.open).abs();
        let prev_body = (prev.close - prev.open).abs();
        if curr_body > prev_body * 1.5 {
            1.0
        } else {
            0.0
        }
    }

    fn compute_gap_indicator(&self) -> f64 {
        if self.bars.len() < 2 {
            return 0.0;
        }
        let curr = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let prev = &self.bars[self.bars.len() - 2];
        if prev.close <= 0.0 || !prev.close.is_finite() {
            return 0.0;
        }
        let gap = curr.open - prev.close;
        safe_clip(gap / prev.close, -0.05, 0.05)
    }

    fn compute_range_position(&self) -> f64 {
        let bar = match self.bars.back() {
            Some(b) => b,
            None => return 0.0,
        };
        let range = bar.high - bar.low + 1e-8;
        safe_clip((bar.close - bar.low) / range, 0.0, 1.0)
    }

    fn compute_volume_momentum(&self, period: usize) -> f64 {
        if self.bars.len() <= period {
            return 0.0;
        }
        let curr_vol = match self.bars.back() {
            Some(b) => b.volume,
            None => return 0.0,
        };
        let prev_vol = self.bars[self.bars.len() - period - 1].volume;
        safe_clip((curr_vol - prev_vol) / (prev_vol + 1e-8), -1.0, 1.0)
    }

    fn compute_volume_acceleration(&self) -> f64 {
        if self.bars.len() < 3 {
            return 0.0;
        }
        let curr = match self.bars.back() {
            Some(b) => b.volume,
            None => return 0.0,
        };
        let prev1 = self.bars[self.bars.len() - 2].volume;
        let prev2 = self.bars[self.bars.len() - 3].volume;
        let vel1 = curr - prev1;
        let vel2 = prev1 - prev2;
        safe_clip(vel1 - vel2, -100.0, 100.0)
    }

    fn compute_volume_max(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        self.bars
            .iter()
            .skip(start)
            .map(|b| b.volume)
            .fold(f64::NEG_INFINITY, f64::max)
    }

    fn compute_volume_min(&self, period: usize) -> f64 {
        let start = self.bars.len().saturating_sub(period);
        self.bars
            .iter()
            .skip(start)
            .map(|b| b.volume)
            .fold(f64::INFINITY, f64::min)
    }

    fn compute_up_down_volume_ratio(&self, period: usize) -> f64 {
        if self.bars.len() < period + 1 {
            return 0.5;
        }
        let start = self.bars.len().saturating_sub(period);
        let mut up_vol = 0.0;
        let mut down_vol = 0.0;
        for i in start..self.bars.len() {
            if i > 0 {
                if self.bars[i].close > self.bars[i - 1].close {
                    up_vol += self.bars[i].volume;
                } else if self.bars[i].close < self.bars[i - 1].close {
                    down_vol += self.bars[i].volume;
                } else {
                    // Unchanged close, volume not attributed to either side
                }
            }
        }
        safe_clip(up_vol / (up_vol + down_vol + 1e-8), 0.0, 1.0)
    }

    fn compute_obv_momentum(&self, period: usize) -> f64 {
        if self.bars.len() < period + 1 {
            return 0.0;
        }
        let mut obv = 0.0;
        let start = self.bars.len().saturating_sub(period);
        for i in (start + 1)..self.bars.len() {
            if self.bars[i].close > self.bars[i - 1].close {
                obv += self.bars[i].volume;
            } else if self.bars[i].close < self.bars[i - 1].close {
                obv -= self.bars[i].volume;
            } else {
                // Unchanged close, OBV unchanged
            }
        }
        safe_clip(obv / 1_000_000.0, -1.0, 1.0)
    }

    fn compute_volume_percentile(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.5;
        }
        let current_vol = match self.bars.back() {
            Some(b) => b.volume,
            None => return 0.5,
        };
        let start = self.bars.len().saturating_sub(period);
        let count_below = self
            .bars
            .iter()
            .skip(start)
            .filter(|b| b.volume < current_vol)
            .count();
        count_below as f64 / period as f64
    }

    fn compute_price_volume_correlation(&self, period: usize) -> f64 {
        if self.bars.len() < period + 1 {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let returns: Vec<f64> = (start + 1..self.bars.len())
            .map(|i| safe_log_return(self.bars[i].close, self.bars[i - 1].close))
            .collect();
        let volumes: Vec<f64> = self.bars.iter().skip(start + 1).map(|b| b.volume).collect();
        self.compute_correlation_from_vecs(&returns, &volumes)
    }

    fn compute_volume_weighted_returns(&self, period: usize) -> f64 {
        if self.bars.len() < period + 1 {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let mut weighted_return = 0.0;
        let mut total_vol = 0.0;
        for i in (start + 1)..self.bars.len() {
            let ret = safe_log_return(self.bars[i].close, self.bars[i - 1].close);
            weighted_return += ret * self.bars[i].volume;
            total_vol += self.bars[i].volume;
        }
        safe_clip(weighted_return / (total_vol + 1e-8), -0.1, 0.1)
    }

    fn compute_range_volume_correlation(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let ranges: Vec<f64> = self
            .bars
            .iter()
            .skip(start)
            .map(|b| (b.high - b.low) / b.close)
            .collect();
        let volumes: Vec<f64> = self.bars.iter().skip(start).map(|b| b.volume).collect();
        self.compute_correlation_from_vecs(&ranges, &volumes)
    }

    fn compute_correlation_from_vecs(&self, x: &[f64], y: &[f64]) -> f64 {
        if x.len() != y.len() || x.is_empty() {
            return 0.0;
        }
        let n = x.len() as f64;
        let mean_x: f64 = x.iter().sum::<f64>() / n;
        let mean_y: f64 = y.iter().sum::<f64>() / n;
        let mut cov = 0.0;
        let mut var_x = 0.0;
        let mut var_y = 0.0;
        for i in 0..x.len() {
            let dx = x[i] - mean_x;
            let dy = y[i] - mean_y;
            cov += dx * dy;
            var_x += dx * dx;
            var_y += dy * dy;
        }
        let denom = (var_x * var_y).sqrt();
        if denom > 1e-8 {
            safe_clip(cov / denom, -1.0, 1.0)
        } else {
            0.0
        }
    }

    fn compute_skewness(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let mean = self.compute_sma(period);
        let std = self.compute_std(period);
        if std < 1e-8 {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let skew: f64 = self
            .bars
            .iter()
            .skip(start)
            .map(|b| ((b.close - mean) / std).powi(3))
            .sum::<f64>()
            / period as f64;
        safe_clip(skew, -3.0, 3.0)
    }

    fn compute_kurtosis(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let mean = self.compute_sma(period);
        let std = self.compute_std(period);
        if std < 1e-8 {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let kurt: f64 = self
            .bars
            .iter()
            .skip(start)
            .map(|b| ((b.close - mean) / std).powi(4))
            .sum::<f64>()
            / period as f64;
        safe_clip(kurt - 3.0, -3.0, 3.0) // Excess kurtosis
    }

    fn compute_percentile(&self, values: &[f64], percentile: f64) -> f64 {
        if values.is_empty() {
            return 0.0;
        }
        let mut sorted = values.to_vec(); // cpu-side sort copy
        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let index = ((sorted.len() as f64 - 1.0) * percentile) as usize;
        sorted[index.min(sorted.len() - 1)]
    }

    fn compute_realized_volatility(&self, period: usize) -> f64 {
        if self.bars.len() < period + 1 {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period + 1);
        let returns: Vec<f64> = (start + 1..self.bars.len())
            .map(|i| safe_log_return(self.bars[i].close, self.bars[i - 1].close))
            .collect();
        let mean = returns.iter().sum::<f64>() / returns.len() as f64;
        let variance =
            returns.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / returns.len() as f64;
        variance.sqrt()
    }

    fn compute_parkinson_volatility(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let sum: f64 = self
            .bars
            .iter()
            .skip(start)
            .map(|b| {
                let hl_ratio = (b.high / b.low).ln();
                hl_ratio * hl_ratio
            })
            .sum();
        (sum / (4.0 * period as f64 * (2.0_f64).ln())).sqrt()
    }

    fn compute_garman_klass_volatility(&self, period: usize) -> f64 {
        if self.bars.len() < period {
            return 0.0;
        }
        let start = self.bars.len().saturating_sub(period);
        let sum: f64 = self
            .bars
            .iter()
            .skip(start)
            .map(|b| {
                let hl = ((b.high / b.low).ln()).powi(2);
                let co = ((b.close / b.open).ln()).powi(2);
                0.5 * hl - (2.0 * (2.0_f64).ln() - 1.0) * co
            })
            .sum();
        (sum / period as f64).sqrt()
    }
}

/// Technical indicator state using common::features shared library
#[derive(Debug)]
struct TechnicalIndicatorState {
    rsi: RSI,
    ema_fast: EMA,
    ema_slow: EMA,
    macd: MACD,
    bollinger: BollingerBands,
    atr: ATR,
    // Cache last computed values
    last_rsi: f64,
    last_ema_fast: f64,
    last_ema_slow: f64,
    last_macd: (f64, f64, f64),
    last_bollinger: (f64, f64, f64),
    last_atr: f64,
}

impl TechnicalIndicatorState {
    fn new() -> Self {
        Self {
            rsi: RSI::new(14),
            ema_fast: EMA::new(12),
            ema_slow: EMA::new(26),
            macd: MACD::new(12, 26, 9),
            bollinger: BollingerBands::new(20, 2.0),
            atr: ATR::new(14),
            last_rsi: 50.0,
            last_ema_fast: 0.0,
            last_ema_slow: 0.0,
            last_macd: (0.0, 0.0, 0.0),
            last_bollinger: (0.0, 0.0, 0.0),
            last_atr: 0.0,
        }
    }

    fn update(&mut self, bar: &OHLCVBar) -> Result<()> {
        // Update all indicators with new bar data and cache results
        self.last_rsi = self.rsi.update(bar.close);
        self.last_ema_fast = self.ema_fast.update(bar.close);
        self.last_ema_slow = self.ema_slow.update(bar.close);
        self.last_macd = self.macd.update(bar.close);
        self.last_bollinger = self.bollinger.update(bar.close);
        self.last_atr = self.atr.update(bar.high, bar.low, bar.close);
        Ok(())
    }
}

// ===== Utility Functions =====

/// Safe log return: log(current / previous), handles edge cases
fn safe_log_return(current: f64, previous: f64) -> f64 {
    if previous <= 0.0 || current <= 0.0 {
        return 0.0;
    }
    let ratio = current / previous;
    if ratio <= 0.0 || !ratio.is_finite() {
        return 0.0;
    }
    ratio.ln()
}

/// Safe normalization: (value - min) / (max - min), clipped to [0, 1]
fn safe_normalize(value: f64, min: f64, max: f64) -> f64 {
    if max <= min || !value.is_finite() {
        return 0.0;
    }
    let normalized = (value - min) / (max - min);
    normalized.clamp(0.0, 1.0)
}

/// Safe clipping: Clip value to [min, max] range
fn safe_clip(value: f64, min: f64, max: f64) -> f64 {
    if !value.is_finite() {
        return 0.0;
    }
    value.clamp(min, max)
}

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

    #[test]
    fn test_feature_extraction_dimensions() {
        // Create synthetic bars
        let bars: Vec<OHLCVBar> = (0..100)
            .map(|i| OHLCVBar {
                timestamp: chrono::Utc::now() + chrono::Duration::hours(i),
                open: 100.0 + i as f64 * 0.1,
                high: 101.0 + i as f64 * 0.1,
                low: 99.0 + i as f64 * 0.1,
                close: 100.5 + i as f64 * 0.1,
                volume: 1000.0 + i as f64 * 10.0,
            })
            .collect();

        let features = extract_ml_features(&bars).unwrap();

        // Should return features for bars after warmup (100 - 50 = 50)
        assert_eq!(features.len(), 50);

        // Each feature vector should be 42-dimensional (40 base + 2 regime: ADX, CUSUM direction)
        for feature_vec in &features {
            assert_eq!(feature_vec.len(), 42);

            // Validate no NaN/Inf
            for &val in feature_vec.iter() {
                assert!(val.is_finite(), "Found non-finite value: {}", val);
            }
        }
    }

    #[test]
    fn test_insufficient_data() {
        let bars: Vec<OHLCVBar> = (0..10)
            .map(|i| OHLCVBar {
                timestamp: chrono::Utc::now() + chrono::Duration::hours(i),
                open: 100.0,
                high: 101.0,
                low: 99.0,
                close: 100.5,
                volume: 1000.0,
            })
            .collect();

        let result = extract_ml_features(&bars);
        assert!(result.is_err());
        assert!(result
            .unwrap_err()
            .to_string()
            .contains("Insufficient data"));
    }

    #[test]
    fn test_safe_log_return() {
        assert_eq!(safe_log_return(110.0, 100.0), (1.1_f64).ln());
        assert_eq!(safe_log_return(0.0, 100.0), 0.0); // Zero current
        assert_eq!(safe_log_return(100.0, 0.0), 0.0); // Zero previous
        assert_eq!(safe_log_return(-10.0, 100.0), 0.0); // Negative
    }

    #[test]
    fn test_safe_normalize() {
        assert_eq!(safe_normalize(50.0, 0.0, 100.0), 0.5);
        assert_eq!(safe_normalize(150.0, 0.0, 100.0), 1.0); // Clipped to 1.0
        assert_eq!(safe_normalize(-50.0, 0.0, 100.0), 0.0); // Clipped to 0.0
        assert_eq!(safe_normalize(f64::NAN, 0.0, 100.0), 0.0); // NaN handling
    }
}
