"""
V29 ORB Structural / Regime Features Diagnostic

Tests three features derived from prior daily bars — orthogonal to OBV-slope / gap-size:

  grav_pull_20_50  : |prior_close - (SMA20+SMA50)/2| / ATR14.
                     High = price far from MA cluster = stretched / strong trend.
                     Low = price near MAs = coiled / stable.
                     Port of pre_event_gravitational_pull (market_features.py) to dict bars.

  market_temp_5_20 : std(5d daily returns) / std(20d daily returns).
                     > 1.0 = vol heating up. < 1.0 = vol cooling (orderly).
                     Port of pre_event_market_temperature (market_features.py).

  momentum_20d     : (prior_close / close_20d_ago) - 1.
                     Positive = rising stock, negative = falling.
                     Distinct from OBV-slope (price return vs volume-weighted trend).

Context: V25 (FINRA short-vol) FAILED. V26 (tape ignition) FAILED.
         V27 (RSI/BB) FAILED — RSI redundant with OBV-slope (ρ=0.726).
         V28 (volatility/compression) FAILED — gap_zscore real but G2 fails (0.18R<0.30R).
         V29: structural positioning axis.

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_grav_pull(prev_bars: list[dict]) -> float | None:
    """Distance from (SMA20+SMA50)/2 normalized by ATR14. Prior day only."""
    if len(prev_bars) < 51:
        return None
    sb = sorted(prev_bars, key=lambda b: b["date"])
    closes = [b["close"] for b in sb]
    sma20 = sum(closes[-20:]) / 20
    sma50 = sum(closes[-50:]) / 50
    ma_center = (sma20 + sma50) / 2
    # ATR14 from last 14 consecutive pairs
    true_ranges = []
    for i in range(max(1, len(sb) - 14), len(sb)):
        curr = sb[i]
        prev = sb[i - 1]
        tr = max(
            curr["high"] - curr["low"],
            abs(curr["high"] - prev["close"]),
            abs(curr["low"] - prev["close"]),
        )
        true_ranges.append(tr)
    if not true_ranges:
        return None
    atr = sum(true_ranges) / len(true_ranges)
    if atr <= 0:
        return None
    return abs(closes[-1] - ma_center) / atr


def compute_market_temp(prev_bars: list[dict]) -> float | None:
    """5d return-std / 20d return-std. < 1.0 = vol cooling."""
    if len(prev_bars) < 21:
        return None
    sb = sorted(prev_bars, key=lambda b: b["date"])[-21:]
    rets = []
    for i in range(1, len(sb)):
        pc = sb[i - 1]["close"]
        cc = sb[i]["close"]
        if pc > 0:
            rets.append((cc - pc) / pc)
    if len(rets) < 20:
        return None

    def _std(xs: list[float]) -> float:
        if len(xs) < 2:
            return 0.0
        m = sum(xs) / len(xs)
        return math.sqrt(sum((x - m) ** 2 for x in xs) / len(xs))

    vol_5d = _std(rets[-5:])
    vol_20d = _std(rets[-20:])
    if vol_20d <= 0:
        return None
    return vol_5d / vol_20d


def compute_momentum_20d(prev_bars: list[dict]) -> float | None:
    """(close_prior / close_21_bars_ago) - 1. Raw 20-day price return."""
    if len(prev_bars) < 21:
        return None
    sb = sorted(prev_bars, key=lambda b: b["date"])
    close_now = sb[-1]["close"]
    close_past = sb[-21]["close"]
    if close_past <= 0:
        return None
    return (close_now - close_past) / close_past


# ── 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
        avg = sum(ys) / len(ys) if ys else 0.0
        return {"n": len(ys), "wr": wr, "avg_r": avg}

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


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

def main() -> None:
    print("=== V29 ORB Structural / Regime 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,
    )
    total_pairs = sum(len(v) for v in candidates.values())
    print(f"Pre-screened: {total_pairs} 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 structural features per trade...")
    sorted_daily: dict[str, list[dict]] = {
        t: sorted(bars, key=lambda b: b["date"]) for t, bars in daily_bars.items()
    }

    annotated: list[dict] = []
    missing: dict[str, int] = {
        "grav_pull_20_50": 0, "market_temp_5_20": 0,
        "momentum_20d": 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_grav = compute_grav_pull(prev_bars)
        f_temp = compute_market_temp(prev_bars)
        f_mom = compute_momentum_20d(prev_bars)
        f_obv = compute_obv_slope_approx(prev_bars, lookback=20)

        enrich_day = enrichment.get(ticker, {}).get(date, {})
        avg_vol = enrich_day.get("avg_daily_vol_14d")

        for fname, fval in [
            ("grav_pull_20_50", f_grav), ("market_temp_5_20", f_temp),
            ("momentum_20d", f_mom), ("obv_slope", f_obv), ("avg_vol", avg_vol),
        ]:
            if fval is None:
                missing[fname] += 1

        annotated.append({
            "ticker": ticker, "date": date, "r": r,
            "grav_pull_20_50": f_grav,
            "market_temp_5_20": f_temp,
            "momentum_20d": f_mom,
            "obv_slope": f_obv,
            "avg_daily_vol": avg_vol,
        })

    total = len(annotated)
    print(f"Annotated: {total} trades")
    for fname in ["grav_pull_20_50", "market_temp_5_20", "momentum_20d"]:
        print(f"  Missing {fname}: {missing[fname]}/{total}")

    feature_defs = [
        ("grav_pull_20_50", "|prior_close − (SMA20+SMA50)/2| / ATR14 — high = stretched from MAs", "low"),
        ("market_temp_5_20", "5d return-std / 20d return-std — < 1 = vol cooling = orderly", "low"),
        ("momentum_20d", "(prior_close / close_20d_ago) − 1 — positive = rising stock", "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(f"  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 V29 config: {weight_sign * weight_magnitude:+.3f}")
    else:
        print(f"\nVERDICT: ABORT — Structural/regime axis null on V24 {LOOKBACK_DAYS}d trade set")
        print("  → V24 remains champion.")
        print("  → Next Ralph iteration: Hurst/OU-theta from libs/features/market_features.py")
        print("    or gap_zscore hard-gate (max_gap_zscore_20d filter) as a filter not a weight")


if __name__ == "__main__":
    main()