"""
speed_study.py -- [2026-09-28] 22-day test of the speed / shape / room measures found
on 9/28 (AGEN, TGB, IMMX vs NWL, SG). Every stock on each day's 9:28 list
(bars/bars_<date>.json.gz), 14-day references as of that morning
(data/reference/<date>/<SYM>.json -- built only from the 14 sessions BEFORE the day).

Measures at every minute (nothing from later in the day):
  move_norm   signed move from the open / the stock's 14-day average move by this minute
  r2_30       smoothness: R^2 of the last 30 closes (signed by slope)
  eff15       straightness: net move / path over 15 bars
  acc15/acc5  acceleration (later half vs earlier half of the steps, in 1-min ATRs)
  vrel15      15-min volume / its 14-day normal for those minutes
  rvol        volume so far / 14-day normal by this minute
  vwap_pct    price vs VWAP
  room_pct    distance to the nearest 14-day resistance band above (None = above all levels)
  pos14       position in the 14-day range (0 = low, 100 = high)

A. Snapshots at 10:00 / 10:30 / 11:00 / 12:00: does each measure separate stocks that
   rise from there to the close (and the day's big winners)? AUC, both halves.
B. Signals (first time per stock-day, 9:45-15:15): what price did next, vs every
   stock at the same minutes (time-matched baseline).
"""
import gzip
import json
import statistics as st
import sys
from datetime import datetime
from pathlib import Path
from zoneinfo import ZoneInfo

HERE = Path(__file__).resolve().parent
sys.path.insert(0, "/var/www/screener/trade")
sys.path.insert(0, "/var/www/screener/trade1/sim/analysis")
import reference as REF
import volatility
from separation import auc

ET = ZoneInfo("America/New_York")


def series(bars, ref):
    t = [datetime.fromtimestamp(b[0], ET) for b in bars]
    idx = [(x.hour - 9) * 60 + x.minute - 30 for x in t]
    o, h, l, c, v = ([b[i] for b in bars] for i in (1, 2, 3, 4, 5))
    day_open = o[0]
    atr, a = [None] * len(bars), None
    pv = vv = 0.0
    vwap = []
    for i in range(len(bars)):
        if i:
            tr = max(h[i] - l[i], abs(h[i] - c[i - 1]), abs(l[i] - c[i - 1]))
            a = tr if a is None else a + (tr - a) / 14
            atr[i] = a
        pv += (h[i] + l[i] + c[i]) / 3 * v[i]
        vv += v[i]
        vwap.append(pv / vv if vv else c[i])
    rows = []
    cum = 0.0
    for i in range(len(bars)):
        cum += v[i]
        m = min(max(idx[i], 0), 389)
        r = {"m": idx[i], "t": t[i].strftime("%H:%M"), "c": c[i],
             "move_norm": (c[i] / day_open - 1) * 100 / max(ref["abs_move_pct"][m], 0.05),
             "rvol": cum / max(ref["cum_vol"][m], 1), "vwap_pct": (c[i] / vwap[i] - 1) * 100}
        for n in (5, 15, 30):
            if i >= n and atr[i]:
                cs = c[i - n:i + 1]
                steps = [cs[k + 1] - cs[k] for k in range(n)]
                net, path = cs[-1] - cs[0], sum(abs(s) for s in steps)
                xm, ym = n / 2, sum(cs) / len(cs)
                sxx = sum((k - xm) ** 2 for k in range(len(cs)))
                slope = sum((k - xm) * (y - ym) for k, y in enumerate(cs)) / sxx
                sst = sum((y - ym) ** 2 for y in cs)
                ssr = sum((y - (ym + slope * (k - xm))) ** 2 for k, y in enumerate(cs))
                half = n // 2
                r[f"r2_{n}"] = (1 - ssr / sst if sst else 0) * (1 if slope > 0 else -1)
                r[f"eff{n}"] = net / path if path else 0
                r[f"acc{n}"] = (sum(steps[half:]) / len(steps[half:]) - sum(steps[:half]) / half) / atr[i]
                norm = sum(ref["vol_per_min"][min(max(x, 0), 389)] for x in idx[i - n + 1:i + 1])
                r[f"vrel{n}"] = sum(v[i - n + 1:i + 1]) / norm if norm else None
        sp = volatility.sr_position(ref.get("sr14") or {}, c[i]) if ref.get("sr14") else {}
        r["room_pct"] = sp.get("dist_to_resistance_pct")
        r["pos14"] = sp.get("pos_in_range_pct")
        r["to_close"] = (c[-1] / c[i] - 1) * 100
        fut = c[i + 1:i + 61]
        r["best60"] = (max(fut) / c[i] - 1) * 100 if fut else None
        r["worst60"] = (min(fut) / c[i] - 1) * 100 if fut else None
        rows.append(r)
    return rows, (c[-1] / day_open - 1) * 100


def load_all():
    data = []
    for p in sorted(HERE.glob("bars/bars_2026-*.json.gz")):
        d = json.load(gzip.open(p, "rt"))
        day = d["date"]
        for s in d["top30_current"]:
            ref = REF.load(day, s)
            bl = d["symbols"][s]["bars"]
            bl = [b for b in bl if datetime.fromtimestamp(b[0], ET).strftime("%H:%M") < "16:00"]
            if ref is None or len(bl) < 100:
                continue
            rows, oc = series(bl, ref)
            data.append({"day": day, "sym": s, "oc": oc, "big": oc >= 5, "rows": rows})
    return data


FEATS = ["move_norm", "r2_30", "r2_15", "eff15", "acc15", "acc5", "vrel15", "rvol", "vwap_pct", "room_pct", "pos14"]


def part_a(data):
    days = sorted({x["day"] for x in data})
    h1 = set(days[:len(days) // 2])
    print(f"\n######## A. Snapshots — {len(data)} stock-days over {len(days)} days "
          f"({sum(x['big'] for x in data)} big winners)")
    for cp in ("10:00", "10:30", "11:00", "12:00"):
        snap = []
        for x in data:
            r = next((r for r in x["rows"] if r["t"] >= cp), None)
            if r and r.get("r2_30") is not None:
                snap.append((x["day"], x["big"], r))
        print(f"\n--- at {cp}: {len(snap)} stock-days; rising to the close {sum(r['to_close'] > 0 for _, _, r in snap)}")
        print(f"   {'measure':10} {'AUC rise':>8} {'1st half':>8} {'2nd half':>8} | {'AUC big':>7} | {'median rose':>11} {'median fell':>11}")
        for f in FEATS:
            up = [r[f] for _, _, r in snap if r["to_close"] > 0 and r.get(f) is not None]
            dn = [r[f] for _, _, r in snap if r["to_close"] <= 0 and r.get(f) is not None]
            u1 = [r[f] for d, _, r in snap if d in h1 and r["to_close"] > 0 and r.get(f) is not None]
            d1 = [r[f] for d, _, r in snap if d in h1 and r["to_close"] <= 0 and r.get(f) is not None]
            u2 = [r[f] for d, _, r in snap if d not in h1 and r["to_close"] > 0 and r.get(f) is not None]
            d2 = [r[f] for d, _, r in snap if d not in h1 and r["to_close"] <= 0 and r.get(f) is not None]
            bg = [r[f] for _, b, r in snap if b and r.get(f) is not None]
            nb = [r[f] for _, b, r in snap if not b and r.get(f) is not None]
            if len(up) > 5 and len(dn) > 5:
                A, A1, A2 = auc(up, dn), auc(u1, d1), auc(u2, d2)
                flag = "  <- consistent" if A1 and A2 and (A1 - .5) * (A2 - .5) > 0 and min(abs(A1 - .5), abs(A2 - .5)) >= .05 else ""
                print(f"   {f:10} {A:8.2f} {A1 or 0:8.2f} {A2 or 0:8.2f} | {auc(bg, nb) or 0:7.2f} | "
                      f"{st.median(up):11.2f} {st.median(dn):11.2f}{flag}")


SIGNALS = {
    "S1 move>=2x + smooth + <=2% over VWAP + room>=1%":
        lambda r: r["move_norm"] >= 2 and (r.get("r2_30") or -1) >= 0.7 and 0 <= r["vwap_pct"] <= 2
        and (r["room_pct"] is None or r["room_pct"] >= 1),
    "S2 move>=2x + smooth": lambda r: r["move_norm"] >= 2 and (r.get("r2_30") or -1) >= 0.7,
    "S3 move>=2x + smooth + <=2% over VWAP": lambda r: r["move_norm"] >= 2 and (r.get("r2_30") or -1) >= 0.7 and 0 <= r["vwap_pct"] <= 2,
    "S4 move>=2x only": lambda r: r["move_norm"] >= 2,
    "S5 smooth only (R2_30>=0.8)": lambda r: (r.get("r2_30") or -1) >= 0.8,
    "S6 AGEN TURN (acc, straight, vol>=2x, <=2% VWAP)":
        lambda r: (r.get("acc15") or -9) >= 0.5 and (r.get("eff15") or 0) >= 0.6 and (r.get("r2_15") or 0) >= 0.7
        and (r.get("vrel15") or 0) >= 2 and 0 <= r["vwap_pct"] <= 2,
    "S7 smooth + into resistance (room < 1%)":
        lambda r: (r.get("r2_30") or -1) >= 0.7 and r["move_norm"] > 0 and r["room_pct"] is not None and r["room_pct"] < 1,
    "CLIMAX (acc5>=1.5, vol>=5x, >=4% VWAP, move>=3x)":
        lambda r: (r.get("acc5") or 0) >= 1.5 and (r.get("vrel15") or 0) >= 5 and r["vwap_pct"] >= 4 and r["move_norm"] >= 3,
}


def part_b(data):
    days = sorted({x["day"] for x in data})
    h1 = set(days[:len(days) // 2])
    base = {}
    for x in data:
        for r in x["rows"]:
            base.setdefault(r["m"], []).append(r["to_close"])
    base = {m: st.mean(v) for m, v in base.items()}
    print(f"\n######## B. Signals (first per stock-day, 9:45-15:15) — to close vs all stocks at the same minute")
    print(f"   {'signal':50} {'events':>6} {'big W':>6} {'rose%':>6} {'avg to close':>12} {'vs base':>8} {'1st/2nd half vs base':>20} {'best60':>7} {'worst60':>8}")
    for name, fn in SIGNALS.items():
        ev = []
        for x in data:
            for r in x["rows"]:
                if 15 <= r["m"] <= 345 and r.get("r2_30") is not None and fn(r):
                    ev.append((x, r))
                    break
        if not ev:
            print(f"   {name:50} {0:6}")
            continue
        ex = [r["to_close"] - base[r["m"]] for _, r in ev]
        e1 = [r["to_close"] - base[r["m"]] for x, r in ev if x["day"] in h1]
        e2 = [r["to_close"] - base[r["m"]] for x, r in ev if x["day"] not in h1]
        print(f"   {name:50} {len(ev):6} {sum(x['big'] for x, _ in ev):6} {sum(r['to_close'] > 0 for _, r in ev) / len(ev) * 100:5.0f}% "
              f"{st.mean(r['to_close'] for _, r in ev):+11.2f}% {st.mean(ex):+7.2f}% "
              f"{(st.mean(e1) if e1 else 0):+9.2f}/{(st.mean(e2) if e2 else 0):+.2f}%  "
              f"{st.mean(r['best60'] for _, r in ev if r['best60'] is not None):+6.2f}% {st.mean(r['worst60'] for _, r in ev if r['worst60'] is not None):+7.2f}%")


if __name__ == "__main__":
    data = load_all()
    part_a(data)
    part_b(data)
    json.dump([{"day": x["day"], "sym": x["sym"], "oc": x["oc"]} for x in data], open(HERE / "speed_study_index.json", "w"))