"""
V30 ORB Trending / Mean-Reversion Features Diagnostic

Tests three features measuring trend persistence and idiosyncratic gap strength:

  hurst_60d       : Hurst exponent via R/S analysis on 60d daily returns.
                    H > 0.5 = trending (persistent), H < 0.5 = mean-reverting.
                    High H before gap-up ORB → stock in trend → better follow-through?
                    Port of pre_event_hurst (market_features.py) to dict bars.

  ou_theta_60d    : OU mean-reversion speed θ = -ln(β) from AR(1) on log prices.
                    Low θ (near 0) = slow reversion = PEAD / breakout friendly.
                    Port of pre_event_ou_theta (market_features.py) to dict bars.

  gap_vs_market   : stock gap_pct − QQQ gap_pct on the same trade day.
                    Measures idiosyncratic (stock-specific) strength above macro open.
                    High = gap driven by stock news/catalyst, not market tide.

Context: V25 short-vol FAILED. V26 tape-ignition FAILED. V27 RSI/BB FAILED (redundant OBV).
         V28 volatility-compression FAILED (gap_zscore G2 near-miss). V29 structural FAILED.
         Note: Hurst+OU previously tested on V23 200d n=96 — near-zero correlation.
               Retesting on V24 400d n=180 for definitive answer.

Gates:
  G1: |Pearson| ≥ 0.07 on ≥ 120 trades (relaxed to 60 if coverage < 80%)
  G2: |top − bottom tercile avg_R| ≥ 0.30R
  G3: |top − bottom tercile WR| ≥ 5pp
  G4: feature coverage ≥ 50%
  G5a: |ρ(feature, obv_slope_20)| < 0.70
  G5b: |ρ(feature, avg_daily_vol_14d)| < 0.70
"""
from __future__ import annotations

import concurrent.futures
import datetime as dt
import math
import os
import sys
from pathlib import Path

import pyarrow.parquet as pq
import yaml
from zoneinfo import ZoneInfo

sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../../..")))

from libs.common.time_utils import trading_days_between
from libs.intraday.domain import ORBStrategyParams
from libs.intraday.features import compute_obv_slope_approx, enrich_daily_bars
from libs.intraday.orb_simulator import ORBSimulationState, run_orb_simulation_with_state
from libs.intraday.screener import orb_pre_screen_candidates

# ── Config ──────────────────────────────────────────────────────────────────
V24_CONFIG = "configs/intraday/strategies/orb_gainers_v24_quality_overlay.yaml"
UNIVERSE_FILE = "configs/symbols_midlarge_snapshot_exact.yaml"
INTRADAY_CACHE_DIR = "data/cache/intraday"
LOOKBACK_DAYS = 400
FEATURE_LOOKBACK_TRADING = 30

_ET = ZoneInfo("America/New_York")
_MKT_OPEN = dt.time(9, 30)
_MKT_CLOSE = dt.time(16, 0)


# ── Feature Computers ────────────────────────────────────────────────────────

def compute_hurst_60d(prev_bars: list[dict]) -> float | None:
    """Hurst exponent via R/S analysis on 60d prior returns."""
    if len(prev_bars) < 62:
        return None
    sb = sorted(prev_bars, key=lambda b: b["date"])[-62:]
    rets = []
    for i in range(1, len(sb)):
        pc, cc = sb[i-1]["close"], sb[i]["close"]
        if pc > 0:
            rets.append((cc - pc) / pc)
    if len(rets) < 30:
        return None

    def rs_stat(series: list[float]) -> float:
        n = len(series)
        mean = sum(series) / n
        devs = [x - mean for x in series]
        cumdev, s = [], 0.0
        for d in devs:
            s += d
            cumdev.append(s)
        r = max(cumdev) - min(cumdev)
        std = (sum(d**2 for d in devs) / n) ** 0.5
        return (r / std) if std > 0 else 0.0

    window_sizes = [w for w in [8, 12, 16, 24, 32] if w <= len(rets) // 2]
    if len(window_sizes) < 2:
        return None

    log_n, log_rs = [], []
    for w in window_sizes:
        rs_vals = []
        for start in range(0, len(rets) - w + 1, w):
            chunk = rets[start:start + w]
            if len(chunk) == w:
                rs_vals.append(rs_stat(chunk))
        if rs_vals:
            avg_rs = sum(rs_vals) / len(rs_vals)
            if avg_rs > 0:
                log_n.append(math.log(w))
                log_rs.append(math.log(avg_rs))

    if len(log_n) < 2:
        return None
    n = len(log_n)
    xm = sum(log_n) / n
    ym = sum(log_rs) / n
    num = sum((log_n[i]-xm)*(log_rs[i]-ym) for i in range(n))
    den = sum((log_n[i]-xm)**2 for i in range(n))
    return (num / den) if den > 0 else 0.5


def compute_ou_theta_60d(prev_bars: list[dict]) -> float | None:
    """OU mean-reversion speed θ = -ln(β) from AR(1) on log prices (60d)."""
    if len(prev_bars) < 61:
        return None
    sb = sorted(prev_bars, key=lambda b: b["date"])[-61:]
    log_prices = []
    for b in sb:
        if b["close"] <= 0:
            return None
        log_prices.append(math.log(b["close"]))
    if len(log_prices) < 61:
        return None
    y = log_prices[1:]
    x = log_prices[:-1]
    n = len(y)
    xm = sum(x) / n
    ym = sum(y) / n
    num = sum((x[i]-xm)*(y[i]-ym) for i in range(n))
    den = sum((x[i]-xm)**2 for i in range(n))
    if den <= 0:
        return None
    beta = num / den
    if beta <= 0 or beta >= 1.0:
        return 0.0
    return -math.log(beta)


# ── Daily Bar Builder ─────────────────────────────────────────────────────────

def _parse_ts(ts_raw: object) -> dt.datetime:
    s = str(ts_raw)
    if s.endswith("Z"):
        s = s[:-1] + "+00:00"
    return dt.datetime.fromisoformat(s).astimezone(_ET)


def _build_daily_bar(path: Path, date: str) -> dict | None:
    try:
        table = pq.read_table(str(path))
        rows = table.to_pydict()
    except Exception:
        return None
    opens, highs, lows, closes, vols = [], [], [], [], []
    for i, ts_raw in enumerate(rows.get("timestamp", [])):
        try:
            ts = _parse_ts(ts_raw)
        except Exception:
            continue
        if _MKT_OPEN <= ts.time() < _MKT_CLOSE:
            opens.append(float(rows["open"][i] or 0))
            highs.append(float(rows["high"][i] or 0))
            lows.append(float(rows["low"][i] or 0))
            closes.append(float(rows["close"][i] or 0))
            vols.append(float(rows["volume"][i] or 0))
    if not opens:
        return None
    return {
        "date": date, "open": opens[0], "high": max(highs),
        "low": min(lows), "close": closes[-1], "volume": sum(vols),
    }


def build_daily_bars(tickers: list[str], dates: list[str], workers: int = 8) -> dict[str, list[dict]]:
    root = Path(INTRADAY_CACHE_DIR)

    def _load(ticker: str) -> tuple[str, list[dict]]:
        d_path = root / ticker
        if not d_path.is_dir():
            return ticker, []
        bars: list[dict] = []
        for date in dates:
            p = d_path / f"{date}.parquet"
            if not p.exists():
                continue
            bar = _build_daily_bar(p, date)
            if bar and bar["close"] > 0:
                bars.append(bar)
        return ticker, bars

    result: dict[str, list[dict]] = {}
    with concurrent.futures.ThreadPoolExecutor(max_workers=workers) as ex:
        for ticker, bars in ex.map(_load, tickers):
            if bars:
                result[ticker] = bars
    return result


def load_intraday_bulk(candidates: dict[str, list[str]]) -> dict[str, dict[str, list[dict]]]:
    import pandas as pd
    result: dict[str, dict[str, list[dict]]] = {}
    for date, tickers in candidates.items():
        day_bars: dict[str, list[dict]] = {}
        for ticker in tickers:
            p = Path(INTRADAY_CACHE_DIR) / ticker / f"{date}.parquet"
            if not p.exists():
                continue
            try:
                df = pd.read_parquet(str(p))
                if not df.empty and len(df) >= 5:
                    day_bars[ticker] = df.to_dict("records")
            except Exception:
                pass
        if day_bars:
            result[date] = day_bars
    return result


# ── Stats Helpers ─────────────────────────────────────────────────────────────

def pearson(xs: list[float], ys: list[float]) -> float | None:
    if len(xs) != len(ys) or len(xs) < 2:
        return None
    n = len(xs)
    xm, ym = sum(xs) / n, sum(ys) / n
    num = sum((xs[i]-xm)*(ys[i]-ym) for i in range(n))
    dx = (sum((x-xm)**2 for x in xs)) ** 0.5
    dy = (sum((y-ym)**2 for y in ys)) ** 0.5
    if dx <= 0 or dy <= 0:
        return None
    return num / (dx * dy)


def tercile_stats(vals: list[float], outcomes_r: list[float]) -> dict:
    if len(vals) < 6:
        return {}
    pairs = sorted(zip(vals, outcomes_r), key=lambda p: p[0])
    n = len(pairs)
    t = n // 3

    def stats(sub):
        ys = [p[1] for p in sub]
        wr = sum(1 for y in ys if y > 0) / len(ys) if ys else 0.0
        return {"n": len(ys), "wr": wr, "avg_r": sum(ys)/len(ys) if ys else 0.0}

    return {"low": stats(pairs[:t]), "mid": stats(pairs[t:2*t]), "high": stats(pairs[2*t:])}


# ── Main ─────────────────────────────────────────────────────────────────────

def main() -> None:
    print("=== V30 ORB Trending / Mean-Reversion Features Diagnostic ===\n")

    with open(V24_CONFIG) as f:
        raw = yaml.safe_load(f)
    params = ORBStrategyParams(**raw["orb_strategy"])
    print(f"V24 config loaded. weight_obv_slope={params.weight_obv_slope}")

    today = dt.date(2026, 4, 21)
    all_td = trading_days_between(today - dt.timedelta(days=700), today)
    trading_days_list = [d.isoformat() for d in all_td[-LOOKBACK_DAYS:]]
    print(f"Window: {trading_days_list[0]} → {trading_days_list[-1]} ({len(trading_days_list)} trading days)")

    extended_td = [d.isoformat() for d in all_td[-(LOOKBACK_DAYS + FEATURE_LOOKBACK_TRADING + 10):]]
    first_cal = dt.date.fromisoformat(extended_td[0])
    last_cal = dt.date.fromisoformat(trading_days_list[-1])
    needed_dates: list[str] = []
    d = first_cal
    while d <= last_cal:
        needed_dates.append(d.isoformat())
        d += dt.timedelta(days=1)

    with open(UNIVERSE_FILE) as f:
        udata = yaml.safe_load(f)
    universe = udata.get("symbols", udata) if isinstance(udata, dict) else udata
    if "QQQ" not in universe:
        universe = list(universe) + ["QQQ"]
    print(f"Universe: {len(universe)} tickers")

    print(f"\nBuilding daily bars ({len(needed_dates)} calendar days)...")
    daily_bars = build_daily_bars(universe, needed_dates)
    print(f"Built daily bars for {len(daily_bars)} tickers")
    print("Computing enrichment...")
    enrichment = enrich_daily_bars(daily_bars, trading_days_list)
    print(f"Enrichment for {len(enrichment)} tickers")

    candidates = orb_pre_screen_candidates(
        daily_bars, trading_days_list, enrichment,
        min_price=params.min_price, min_atr=params.min_atr_14,
        min_avg_dollar_vol=params.min_avg_dollar_volume, max_per_day=None,
    )
    print(f"Pre-screened: {sum(len(v) for v in candidates.values())} ticker-days")
    print("Loading intraday bars...")
    all_intraday = load_intraday_bulk(candidates)
    print(f"Loaded: {sum(len(v) for v in all_intraday.values())} ticker-days")

    print("\nRunning V24 simulation...")
    state = ORBSimulationState(equity=params.initial_capital)
    day_results, _ = run_orb_simulation_with_state(
        all_intraday, trading_days_list, params, enrichment, state=state,
    )
    all_trades = [t for dr in day_results for t in dr.trades]
    trades_with_r = [t for t in all_trades if getattr(t, "r_multiple_at_exit", None) is not None]
    print(f"Total trades: {len(all_trades)}, with r_multiple: {len(trades_with_r)}")

    if len(trades_with_r) < 20:
        print("ABORT: fewer than 20 trades with r_multiple")
        return

    print("\nComputing trending/mean-reversion features per trade...")
    sorted_daily: dict[str, list[dict]] = {
        t: sorted(bars, key=lambda b: b["date"]) for t, bars in daily_bars.items()
    }

    # Pre-build QQQ gap_pct lookup per trade day from enrichment
    qqq_enrich = enrichment.get("QQQ", {})

    annotated: list[dict] = []
    missing: dict[str, int] = {
        "hurst_60d": 0, "ou_theta_60d": 0, "gap_vs_market": 0,
        "obv_slope": 0, "avg_vol": 0,
    }

    for trade in trades_with_r:
        ticker = trade.ticker
        date = str(trade.date)[:10]
        r = float(trade.r_multiple_at_exit)

        bars_t = sorted_daily.get(ticker, [])
        prev_bars = [b for b in bars_t if b["date"][:10] < date]

        f_hurst = compute_hurst_60d(prev_bars)
        f_ou = compute_ou_theta_60d(prev_bars)
        f_obv = compute_obv_slope_approx(prev_bars, lookback=20)

        # gap_vs_market: stock gap_pct - QQQ gap_pct on trade day
        enrich_day = enrichment.get(ticker, {}).get(date, {})
        stock_gap = enrich_day.get("gap_pct")
        qqq_day_enrich = qqq_enrich.get(date, {})
        qqq_gap = qqq_day_enrich.get("gap_pct")
        f_gap_vs_mkt = (stock_gap - qqq_gap) if (stock_gap is not None and qqq_gap is not None) else None

        avg_vol = enrich_day.get("avg_daily_vol_14d")

        for fname, fval in [
            ("hurst_60d", f_hurst), ("ou_theta_60d", f_ou),
            ("gap_vs_market", f_gap_vs_mkt), ("obv_slope", f_obv), ("avg_vol", avg_vol),
        ]:
            if fval is None:
                missing[fname] += 1

        annotated.append({
            "ticker": ticker, "date": date, "r": r,
            "hurst_60d": f_hurst,
            "ou_theta_60d": f_ou,
            "gap_vs_market": f_gap_vs_mkt,
            "obv_slope": f_obv,
            "avg_daily_vol": avg_vol,
        })

    total = len(annotated)
    print(f"Annotated: {total} trades")
    for fname in ["hurst_60d", "ou_theta_60d", "gap_vs_market"]:
        print(f"  Missing {fname}: {missing[fname]}/{total}")

    feature_defs = [
        ("hurst_60d", "R/S Hurst exponent on 60d returns — H>0.5 = trending", "high"),
        ("ou_theta_60d", "OU θ = -ln(β) from AR(1) — low = slow reversion = breakout-friendly", "low"),
        ("gap_vs_market", "stock gap_pct − QQQ gap_pct — high = idiosyncratic strength", "high"),
    ]

    print("\n" + "=" * 90)
    print(f"FEATURE ANALYSIS — V24 {LOOKBACK_DAYS}d trade set")
    print("=" * 90)

    results: dict[str, dict | None] = {}
    for feat_name, description, hypothesized_best in feature_defs:
        valid = [(a[feat_name], a["r"]) for a in annotated if a[feat_name] is not None]
        if len(valid) < 20:
            print(f"\n{feat_name}: SKIP — only {len(valid)} valid trades (need ≥20)")
            results[feat_name] = None
            continue

        vals = [v[0] for v in valid]
        rs = [v[1] for v in valid]
        overall_wr = sum(1 for v in valid if v[1] > 0) / len(valid)
        rho = pearson(vals, rs)
        tstat = tercile_stats(vals, rs)
        coverage_pct = len(valid) / total
        rho_abs = abs(rho) if rho is not None else 0.0

        if tstat:
            best_tercile = "high" if tstat["high"]["avg_r"] >= tstat["low"]["avg_r"] else "low"
            worst_tercile = "low" if best_tercile == "high" else "high"
            direction_match = best_tercile == hypothesized_best

            print(f"\n{'─'*60}")
            print(f"FEATURE: {feat_name}")
            print(f"  Description: {description}")
            print(f"  n={len(valid)}/{total} ({coverage_pct*100:.0f}% coverage), overall WR={overall_wr*100:.1f}%")
            print(f"  Pearson(feature, r_multiple) = {rho:.4f}" if rho is not None else "  Pearson = n/a")
            print("  Tercile breakdown (low→high feature value):")
            h, m, lo = tstat["high"], tstat["mid"], tstat["low"]
            print(f"    Bottom: n={lo['n']}, WR={lo['wr']*100:.1f}%, avg_R={lo['avg_r']:+.3f}")
            print(f"    Middle: n={m['n']}, WR={m['wr']*100:.1f}%, avg_R={m['avg_r']:+.3f}")
            print(f"    Top:    n={h['n']}, WR={h['wr']*100:.1f}%, avg_R={h['avg_r']:+.3f}")
            print(f"  Empirical best: '{best_tercile}' tercile (hypothesis: '{hypothesized_best}' → {'✓ confirmed' if direction_match else '✗ INVERTED'})")

            best = tstat[best_tercile]
            worst = tstat[worst_tercile]
            min_trades = 120 if coverage_pct >= 0.80 else 60
            g1 = rho_abs >= 0.07 and len(valid) >= min_trades
            g2 = best["avg_r"] - worst["avg_r"] >= 0.30
            g3 = best["wr"] >= worst["wr"] + 0.05
            g4 = coverage_pct >= 0.50
            n_failed = sum([not g1, not g2, not g3, not g4])
            overall_pass = n_failed == 0

            print(f"  Gates (best='{best_tercile}' tercile):")
            print(f"    G1 |Pearson|≥0.07 + n≥{min_trades}: {rho_abs:.4f}, n={len(valid)}  → {'PASS ✓' if g1 else 'FAIL ✗'}")
            print(f"    G2 avg_R gap ≥ 0.30R: {best['avg_r'] - worst['avg_r']:+.3f}  → {'PASS ✓' if g2 else 'FAIL ✗'}")
            print(f"    G3 WR gap ≥ 5pp: {(best['wr'] - worst['wr'])*100:+.1f}pp  → {'PASS ✓' if g3 else 'FAIL ✗'}")
            print(f"    G4 coverage ≥ 50%: {coverage_pct*100:.0f}%  → {'PASS ✓' if g4 else 'FAIL ✗'}")
            print(f"  VERDICT: {'ALL GATES PASS → PROCEED TO PHASE 2' if overall_pass else f'FAIL ({n_failed} gate(s) failed)'}")

            results[feat_name] = {
                "pass": overall_pass, "pearson": rho, "stats": tstat,
                "n": len(valid), "coverage": coverage_pct, "best_tercile": best_tercile,
            }

    print(f"\n{'─'*60}")
    print("PAIRWISE CORRELATIONS:")
    for feat_name in [fd[0] for fd in feature_defs]:
        for ref_key, ref_label, gate_name in [
            ("obv_slope", "obv_slope_20", "G5a"),
            ("avg_daily_vol", "avg_daily_vol_14d", "G5b"),
        ]:
            combined = [
                (a[feat_name], a[ref_key]) for a in annotated
                if a[feat_name] is not None and a[ref_key] is not None
            ]
            if len(combined) >= 10:
                rho_x = pearson([c[0] for c in combined], [c[1] for c in combined])
                if rho_x is not None:
                    gate_pass = abs(rho_x) < 0.70
                    print(f"  ρ({feat_name[:28]:28s}, {ref_label}): {rho_x:+.4f}  {gate_name}: {'PASS ✓' if gate_pass else 'FAIL ✗'}")

    print("\n  Inter-feature correlations:")
    feat_keys = [fd[0] for fd in feature_defs]
    for i, f1 in enumerate(feat_keys):
        for f2 in feat_keys[i+1:]:
            combined = [(a[f1], a[f2]) for a in annotated if a[f1] is not None and a[f2] is not None]
            if len(combined) >= 10:
                rho_x = pearson([c[0] for c in combined], [c[1] for c in combined])
                if rho_x is not None:
                    print(f"    ρ({f1[:28]:28s}, {f2[:28]:28s}): {rho_x:+.4f}")

    passing = [name for name, r in results.items() if r is not None and r["pass"]]
    failed = [name for name, r in results.items() if r is not None and not r["pass"]]
    skipped = [name for name, r in results.items() if r is None]

    print(f"\n{'='*90}")
    print("SUMMARY")
    print(f"{'='*90}")
    print(f"Features passing all gates: {passing if passing else 'NONE'}")
    print(f"Features failing gates:     {failed if failed else 'NONE'}")
    print(f"Features skipped:           {skipped if skipped else 'NONE'}")

    if passing:
        best = max(passing, key=lambda n: abs(results[n]["pearson"] or 0))
        pearson_val = results[best]["pearson"]
        obv_pearson_ref = 0.2349
        weight_magnitude = round(0.05 * min(1.0, abs(pearson_val) / obv_pearson_ref), 2)
        weight_magnitude = max(weight_magnitude, 0.02)
        best_direction = results[best]["best_tercile"]
        weight_sign = +1 if best_direction == "high" else -1

        print(f"\nVERDICT: PROCEED TO PHASE 2")
        print(f"  Best feature: {best}")
        print(f"  Pearson: {pearson_val:.4f}")
        print(f"  Direction: '{best_direction}' tercile is best → weight = {weight_sign * weight_magnitude:+.3f}")
        print(f"  Suggested weight in V30 config: {weight_sign * weight_magnitude:+.3f}")
    else:
        print(f"\nVERDICT: ABORT — Trending/mean-reversion axis null on V24 {LOOKBACK_DAYS}d trade set")
        print("  → V24 remains champion.")
        print("  → Signal inventory exhausted from libs/features/market_features.py.")
        print("  → Next Ralph iteration: gap_zscore hard-gate backtest (max_gap_zscore_20d filter)")
        print("    or momentum_20d + grav_pull composite weight (both near-miss, orthogonal)")


if __name__ == "__main__":
    main()