"""What do the system's OWN indicators read right before a violent bar?  Analysis tooling only.

Uses backtest.gpu_backtest.precompute_indicators (the same series the backtester feeds the
engine) on 3 years of bars; every state is taken at t-1 (known before the bar). Reports burst
rate (range >= 3x prior ATR14) for the NEXT bar and within the NEXT HOUR by indicator state,
by regime, by volatility regime, by weekday x hour, day of month, calendar proximity — each in
both halves — and whether the state predicts the burst's DIRECTION.

    python backtest/research/burst_indicator_map.py PARQUET BAR_MINUTES [--start --end]
"""
import argparse, sys
from pathlib import Path
import numpy as np, pandas as pd
ROOT = Path(__file__).resolve().parents[2]; sys.path.insert(0, str(ROOT)); sys.path.insert(0, str(ROOT / "backtest"))
from gpu_backtest import precompute_indicators
from config import TechnicalConfig

CAL = ROOT / "backtest" / "cache" / "calendar.parquet"


def wilder_atr(df, period=14):
    pc = df["Close"].shift(1)
    tr = pd.concat([df["High"] - df["Low"], (df["High"] - pc).abs(), (df["Low"] - pc).abs()], axis=1).max(axis=1)
    return tr.ewm(alpha=1 / period, min_periods=period, adjust=False).mean()


def rate_table(d, col, target, half, order=None):
    base1 = d.loc[half == 1, target].mean(); base2 = d.loc[half == 2, target].mean()
    g = d.groupby(col, observed=True)
    out = pd.DataFrame({"n": g.size(),
                        "h1_lift": d[half == 1].groupby(col, observed=True)[target].mean() / base1,
                        "h2_lift": d[half == 2].groupby(col, observed=True)[target].mean() / base2,
                        "rate%": 100 * g[target].mean()})
    if order is not None: out = out.reindex([o for o in order if o in out.index])
    out["stable"] = np.where((out.h1_lift > 1.15) & (out.h2_lift > 1.15), "HIGH", np.where((out.h1_lift < 0.85) & (out.h2_lift < 0.85), "low", ""))
    return out.round(2)


def main():
    ap = argparse.ArgumentParser(); ap.add_argument("parquet"); ap.add_argument("bar_minutes", type=int)
    ap.add_argument("--start", default="2023-09-05"); ap.add_argument("--end", default="2026-09-05")
    a = ap.parse_args(); bm = a.bar_minutes
    px = pd.read_parquet(a.parquet)
    px = px[(px.index >= pd.Timestamp(a.start, tz="UTC") - pd.Timedelta(days=30)) & (px.index < pd.Timestamp(a.end, tz="UTC"))]
    ic = precompute_indicators(px, TechnicalConfig())
    d = pd.DataFrame(index=px.index)
    atr_prev = wilder_atr(px).shift(1)
    rng = px["High"] - px["Low"]
    d["burst"] = (rng >= 3.0 * atr_prev).astype(int)
    H = max(1, 60 // bm)
    d["next_hour"] = d["burst"][::-1].rolling(H, min_periods=1).max()[::-1].astype(int)
    d["next_bar"] = d["burst"]
    d["bdir"] = np.sign(px["Close"] - px["Open"])
    # states at t-1
    sh = lambda arr: pd.Series(np.asarray(arr, dtype=float), index=px.index).shift(1)
    rsi = sh(ic.rsi); adx = sh(ic.adx); adxs = sh(ic.adx_slope); hist = sh(ic.macd_hist); bbp = sh(ic.bb_position)
    bbw = sh(ic.bb_band_width); bbaw = sh(ic.bb_avg_width); kcsq = sh(ic.kc_squeeze); roc = sh(ic.roc)
    ef, em, es = sh(ic.ema_fast), sh(ic.ema_medium), sh(ic.ema_slow); et = sh(ic.ema_trend) if ic.ema_trend is not None else None
    close1 = px["Close"].shift(1); atr1 = sh(ic.atr)
    kcu, kcl = sh(ic.kc_upper), sh(ic.kc_lower)
    d["rsi"] = pd.cut(rsi, [0, 30, 40, 50, 60, 70, 100], labels=["<30", "30-40", "40-50", "50-60", "60-70", ">=70"])
    d["adx"] = pd.cut(adx, [0, 20, 25, 35, 50, 200], labels=["<20", "20-25", "25-35", "35-50", ">=50"])
    d["adx_slope"] = pd.cut(adxs, [-1e9, -3, 0, 3, 1e9], labels=["falling<-3", "-3..0", "0..3", "rising>3"])
    hist_prev = hist.shift(1)
    d["macd"] = np.select([(hist > 0) & (hist > hist_prev), (hist > 0), (hist < 0) & (hist < hist_prev), (hist < 0)],
                          ["pos&rising", "pos&fading", "neg&falling", "neg&fading"], default="flat")
    d["bb_pos"] = pd.cut(bbp, [-1e9, 0, 0.2, 0.8, 1.0, 1e9], labels=["<0 (below)", "0-0.2", "0.2-0.8", "0.8-1", ">1 (above)"])
    d["bb_squeeze"] = np.where(bbw < 0.75 * bbaw, "squeeze", "normal")
    d["kc_squeeze"] = np.where(kcsq > 0, "TTM squeeze", "no")
    d["kc_pos"] = np.select([close1 > kcu, close1 < kcl], ["above KC", "below KC"], default="inside")
    align = np.sign(close1 - ef) + np.sign(ef - em) + np.sign(em - es) + (np.sign(es - et) if et is not None else 0)
    d["ema_align"] = pd.cut(align, [-5, -3.5, -0.5, 0.5, 3.5, 5], labels=["fully bear(-4)", "mostly bear", "mixed", "mostly bull", "fully bull(+4)"])
    ribbon = (pd.concat([ef, em, es], axis=1).max(axis=1) - pd.concat([ef, em, es], axis=1).min(axis=1)) / atr1
    d["ribbon"] = np.where(ribbon < 1.0, "compressed(<1 ATR)", "expanded")
    ema_dir = np.where(align > 0, 1, np.where(align < 0, -1, 0))
    d["regime"] = np.select([(adx >= 25) & (ema_dir > 0), (adx >= 25) & (ema_dir < 0), adx >= 25, adxs > 3.0],
                            ["trending_up", "trending_down", "trending(mixed)", "transitioning"], default="ranging")
    atr_pct = atr1.rolling(50).apply(lambda w: (w[-1] > w).mean() * 100, raw=True)
    d["vol_regime"] = pd.cut(atr_pct, [-1, 30, 60, 80, 101], labels=["low", "normal", "high", "extreme"])
    d["roc"] = pd.cut(roc, [-1e9, -1, -0.3, 0.3, 1, 1e9], labels=["<-1%", "-1..-0.3", "-0.3..0.3", "0.3..1", ">1%"])
    d["ema200_dist"] = pd.cut((close1 - et) / et * 100 if et is not None else pd.Series(np.nan, index=px.index),
                              [-1e9, -2.5, -0.5, 0.5, 2.5, 1e9], labels=["<-2.5%", "-2.5..-0.5", "near(±0.5%)", "0.5..2.5", ">2.5%"])
    # time
    ct = px.index + pd.Timedelta(minutes=bm)
    d["hour"] = ct.hour; d["weekday"] = ct.weekday; d["dom"] = pd.cut(ct.day, [0, 5, 10, 15, 20, 25, 31], labels=["1-5", "6-10", "11-15", "16-20", "21-25", "26-31"])
    d["first_friday"] = np.where((ct.weekday == 4) & (ct.day <= 7), "first Friday", "other")
    # calendar proximity (high-impact USD)
    cal = pd.read_parquet(CAL); cal["ts"] = pd.to_datetime(cal.dateline, unit="s", utc=True)
    hi = np.sort(cal[(cal.currency == "USD") & (cal.impactName.str.lower() == "high")].ts.values.astype("datetime64[ns]"))
    t = ct.values.astype("datetime64[ns]")
    nxt = np.searchsorted(hi, t)  # next event at/after bar close
    mins_to = np.where(nxt < len(hi), (hi[np.minimum(nxt, len(hi) - 1)] - t) / np.timedelta64(1, "m"), 1e9)
    prv = nxt - 1
    mins_since = np.where(prv >= 0, (t - hi[np.maximum(prv, 0)]) / np.timedelta64(1, "m"), 1e9)
    d["calendar"] = np.select([mins_since <= 30, mins_to <= 15, mins_to <= 60, mins_to <= 180], ["0-30m AFTER hi-impact", "next 15m", "15-60m ahead", "1-3h ahead"], default="none within 3h")
    d = d[d.index >= pd.Timestamp(a.start, tz="UTC")].dropna(subset=["rsi", "adx", "bb_pos", "vol_regime"])
    half = pd.Series(np.where(np.arange(len(d)) < len(d) // 2, 1, 2), index=d.index)
    print(f"bars {len(d)}  base next_bar {100*d.next_bar.mean():.3f}%  next_hour {100*d.next_hour.mean():.2f}%  (lift = rate / base; 'HIGH' = >1.15x in BOTH halves, 'low' = <0.85x in both)")
    groups = ["calendar", "regime", "vol_regime", "rsi", "adx", "adx_slope", "macd", "bb_pos", "bb_squeeze", "kc_squeeze", "kc_pos", "ema_align", "ribbon", "roc", "ema200_dist", "dom", "first_friday"]
    for target in ("next_bar", "next_hour"):
        print(f"\n################ TARGET: {target} ################")
        for g in groups:
            print(f"-- {g}"); print(rate_table(d, g, target, half).to_string())
    # weekday x hour top cells (next_bar)
    piv = d.pivot_table(index="weekday", columns="hour", values="next_bar", aggfunc="mean") * 100
    base = 100 * d.next_bar.mean()
    cells = piv.stack().sort_values(ascending=False)
    print(f"\n################ TOP 15 weekday x close-hour cells (next-bar burst pct, base {base:.2f}) ################")
    names = {0: "Mon", 1: "Tue", 2: "Wed", 3: "Thu", 4: "Fri", 6: "Sun"}
    for (wd, h), v in cells.head(15).items():
        n = int(((d.weekday == wd) & (d.hour == h)).sum())
        print(f"  {names.get(wd, wd)} {h:02d}:xx  {v:.2f}%  ({v/base:.1f}x, n={n})")
    # direction: does the state predict the burst direction?
    b = d[d.burst == 1]
    print("\n################ BURST DIRECTION vs prior state (share of bursts going UP) ################")
    for g, up_states in (("ema_align", None), ("rsi", None), ("macd", None), ("regime", None)):
        print(f"-- {g}: " + ", ".join(f"{k}: {100*(v>0).mean():.0f}% up (n={len(v)})" for k, v in b.groupby(g, observed=True).bdir if len(v) >= 20))


if __name__ == "__main__":
    main()