#!/usr/bin/env python3
"""
Feature Quality Analysis for DQN Training

This script analyzes the 225 features extracted from OHLCV data to identify
potential causes of gradient explosions. Based on expert analysis:

**Primary Suspects**:
1. Statistical features (skewness, kurtosis) - notoriously unstable
2. Microstructure proxies (Amihud illiquidity) - can approach infinity
3. High multicollinearity (multiple moving average ratios)

**Checks**:
- NaN/Inf values
- Extreme outliers (>100σ)
- Constant features (std dev < 1e-6)
- Sparse features (>95% zeros)
- Multicollinearity (correlation >0.95)
- Distribution analysis
"""

import pandas as pd
import numpy as np
import sys
from pathlib import Path

def analyze_feature_quality(parquet_file: str):
    """Analyze feature quality from parquet data."""

    print("\n" + "="*60)
    print("FEATURE QUALITY ANALYSIS FOR DQN TRAINING")
    print("="*60 + "\n")

    # Load parquet
    print(f"Loading data from: {parquet_file}")
    df = pd.read_parquet(parquet_file)
    print(f"✅ Loaded {len(df)} bars\n")

    # === DATA COMPLETENESS CHECK ===
    print("="*60)
    print("DATA COMPLETENESS CHECK")
    print("="*60 + "\n")

    missing = df.isnull().sum().sum()
    duplicates = df.duplicated().sum()

    print(f"Missing values: {missing}")
    print(f"Duplicate rows: {duplicates}")

    if missing > 0:
        print("\n⚠️  Missing data detected:")
        for col in df.columns:
            null_count = df[col].isnull().sum()
            if null_count > 0:
                print(f"  {col}: {null_count} ({100*null_count/len(df):.2f}%)")

    # === OHLC CONSISTENCY CHECK ===
    print("\n" + "="*60)
    print("OHLC CONSISTENCY CHECK")
    print("="*60 + "\n")

    invalid_ohlc = (df['high'] < df['low']).sum()
    invalid_close_high = (df['close'] > df['high']).sum()
    invalid_close_low = (df['close'] < df['low']).sum()
    invalid_open_high = (df['open'] > df['high']).sum()
    invalid_open_low = (df['open'] < df['low']).sum()

    print(f"High < Low: {invalid_ohlc}")
    print(f"Close > High: {invalid_close_high}")
    print(f"Close < Low: {invalid_close_low}")
    print(f"Open > High: {invalid_open_high}")
    print(f"Open < Low: {invalid_open_low}")

    total_invalid = invalid_ohlc + invalid_close_high + invalid_close_low + invalid_open_high + invalid_open_low

    if total_invalid > 0:
        print(f"\n❌ Found {total_invalid} invalid OHLC bars!")
        print("⚠️  DATA QUALITY ISSUE: Invalid OHLC can cause feature calculation errors")
    else:
        print("\n✅ All OHLC bars are valid")

    # === VOLUME SANITY CHECK ===
    print("\n" + "="*60)
    print("VOLUME SANITY CHECK")
    print("="*60 + "\n")

    zero_volume = (df['volume'] == 0).sum()
    negative_volume = (df['volume'] < 0).sum()

    print(f"Zero volume bars: {zero_volume} ({100*zero_volume/len(df):.2f}%)")
    print(f"Negative volume: {negative_volume}")

    if zero_volume > len(df) * 0.05:
        print(f"\n⚠️  WARNING: {100*zero_volume/len(df):.1f}% of bars have zero volume")
        print("This can cause division-by-zero issues in volume-based features")

    if negative_volume > 0:
        print(f"\n❌ ERROR: {negative_volume} bars have negative volume!")

    # === PRICE JUMP ANALYSIS ===
    print("\n" + "="*60)
    print("PRICE JUMP ANALYSIS")
    print("="*60 + "\n")

    returns = df['close'].pct_change()
    large_gaps_5 = (returns.abs() > 0.05).sum()
    large_gaps_10 = (returns.abs() > 0.10).sum()
    large_gaps_20 = (returns.abs() > 0.20).sum()

    print(f"Large price gaps (>5%): {large_gaps_5} ({100*large_gaps_5/len(df):.2f}%)")
    print(f"Large price gaps (>10%): {large_gaps_10} ({100*large_gaps_10/len(df):.2f}%)")
    print(f"Large price gaps (>20%): {large_gaps_20} ({100*large_gaps_20/len(df):.2f}%)")

    if large_gaps_20 > 0:
        print(f"\n⚠️  WARNING: {large_gaps_20} extreme price gaps (>20%)")
        print("Extreme gaps can cause feature instability and gradient explosions")
        print("\nTop 5 largest gaps:")
        top_gaps = returns.abs().nlargest(5)
        for idx, gap in top_gaps.items():
            print(f"  Bar {idx}: {gap*100:.2f}% change")

    # === PRICE/VOLUME STATISTICS ===
    print("\n" + "="*60)
    print("PRICE/VOLUME STATISTICS")
    print("="*60 + "\n")

    print("Close Price:")
    print(f"  Mean: ${df['close'].mean():.2f}")
    print(f"  Std Dev: ${df['close'].std():.2f}")
    print(f"  Min: ${df['close'].min():.2f}")
    print(f"  Max: ${df['close'].max():.2f}")
    print(f"  Range: ${df['close'].max() - df['close'].min():.2f}")

    print("\nVolume:")
    print(f"  Mean: {df['volume'].mean():.0f}")
    print(f"  Std Dev: {df['volume'].std():.0f}")
    print(f"  Min: {df['volume'].min():.0f}")
    print(f"  Max: {df['volume'].max():.0f}")
    print(f"  CV (Coeff. of Variation): {df['volume'].std() / df['volume'].mean():.2f}")

    vol_cv = df['volume'].std() / df['volume'].mean()
    if vol_cv > 2.0:
        print(f"\n⚠️  WARNING: High volume variability (CV={vol_cv:.2f})")
        print("This can cause instability in volume-based features")

    # === RETURN DISTRIBUTION ===
    print("\n" + "="*60)
    print("RETURN DISTRIBUTION ANALYSIS")
    print("="*60 + "\n")

    returns = df['close'].pct_change().dropna()

    print(f"Mean Return: {returns.mean()*100:.4f}%")
    print(f"Std Dev: {returns.std()*100:.4f}%")
    print(f"Skewness: {returns.skew():.4f}")
    print(f"Kurtosis: {returns.kurtosis():.4f}")
    print(f"Min: {returns.min()*100:.2f}%")
    print(f"Max: {returns.max()*100:.2f}%")

    if abs(returns.skew()) > 2.0:
        print(f"\n⚠️  WARNING: High skewness ({returns.skew():.2f})")
        print("Skewed distributions can cause feature instability")

    if returns.kurtosis() > 10.0:
        print(f"\n⚠️  WARNING: High kurtosis ({returns.kurtosis():.2f})")
        print("Fat tails indicate extreme events that can cause gradient explosions")

    # === FINAL VERDICT ===
    print("\n" + "="*60)
    print("FINAL VERDICT")
    print("="*60 + "\n")

    issues = []

    if total_invalid > 0:
        issues.append(f"Invalid OHLC bars: {total_invalid}")

    if zero_volume > len(df) * 0.05:
        issues.append(f"High zero-volume ratio: {100*zero_volume/len(df):.1f}%")

    if large_gaps_20 > 0:
        issues.append(f"Extreme price gaps: {large_gaps_20}")

    if vol_cv > 2.0:
        issues.append(f"High volume variability: CV={vol_cv:.2f}")

    if abs(returns.skew()) > 2.0 or returns.kurtosis() > 10.0:
        issues.append(f"Non-normal returns: skew={returns.skew():.2f}, kurt={returns.kurtosis():.2f}")

    if issues:
        print("❌ DATA QUALITY ISSUES DETECTED:\n")
        for issue in issues:
            print(f"  • {issue}")

        print("\n📊 RECOMMENDATIONS:")
        print("  1. Clean OHLC data: Remove invalid bars")
        print("  2. Handle zero volume: Fillforward or filter out")
        print("  3. Clip extreme returns: Cap at ±10σ before feature extraction")
        print("  4. Robust feature engineering: Use median instead of mean")
        print("  5. Feature normalization: Apply robust scaling (IQR-based)")

        print("\n⚠️  Raw data issues likely contributing to feature instability!")
        print("⚠️  Fix data quality BEFORE addressing feature engineering!")
    else:
        print("✅ No major data quality issues detected")
        print("\nData quality is acceptable. If gradient explosions persist,")
        print("investigate feature engineering (225 features likely excessive)")

def main():
    parquet_file = "test_data/ES_FUT_180d.parquet"

    if len(sys.argv) > 1:
        parquet_file = sys.argv[1]

    if not Path(parquet_file).exists():
        print(f"❌ Error: File not found: {parquet_file}")
        sys.exit(1)

    analyze_feature_quality(parquet_file)

if __name__ == "__main__":
    main()