#!/usr/bin/env python3
"""
shapelab.py -- "Why can't we fit math to the SHAPE, noise or not?"

The radiology protocol. A cancer detector isn't validated by "does it find patterns"
-- it's validated by whether it DISCRIMINATES against ground truth, out of sample,
versus a control. So: train a real gradient-boosted classifier on multi-bar shape
features, measure out-of-sample AUC, and run the identical pipeline on shuffled
(random-walk) prices as the control arm.

  AUC_real ~= AUC_shuffled ~= 0.50  -> the shape carries no information
  AUC_real >  AUC_shuffled          -> there IS learnable structure

Run against TWO targets, because they are not the same question:
  DIRECTION  -- which way does price go next?
  VOLATILITY -- how MUCH does it move next? (the "flat parts vs active parts")
"""
import glob, os
import numpy as np
import pandas as pd
from sklearn.metrics import roc_auc_score
import lightgbm as lgb

FUT = "user_data/data/krakenfutures/futures"
SEED = 20260728
HORIZON = 4          # predict 4h ahead


def features(c, v):
    """Causal shape descriptors -- everything uses data up to and including bar t."""
    c = pd.Series(c); v = pd.Series(v)
    lr = np.log(c).diff()
    X = pd.DataFrame(index=c.index)
    for k in (1, 2, 4, 8, 24, 72, 168):
        X[f"ret{k}"] = np.log(c).diff(k)
    for k in (6, 24, 72, 168):
        X[f"vol{k}"] = lr.rolling(k).std()
        X[f"absret{k}"] = lr.abs().rolling(k).mean()
    for k in (24, 72, 168):
        hi = c.rolling(k).max(); lo = c.rolling(k).min()
        X[f"pos{k}"] = (c - lo) / (hi - lo)                 # where in the range
        X[f"rng{k}"] = (hi - lo) / c                        # how wide the range
    for k in (12, 48):
        m = c.rolling(k).mean()
        X[f"dev{k}"] = (c - m) / m                          # extension from mean
        X[f"slope{k}"] = m.diff(k) / m                      # trend of the mean
    # momentum-of-vol, vol-of-vol: does the regime itself have shape?
    X["volratio"] = lr.rolling(6).std() / lr.rolling(72).std()
    X["volvol"] = lr.rolling(24).std().rolling(72).std()
    X["skew72"] = lr.rolling(72).skew()
    X["kurt72"] = lr.rolling(72).kurt()
    X["vz24"] = (v - v.rolling(24).mean()) / v.rolling(24).std()
    X["vz168"] = (v - v.rolling(168).mean()) / v.rolling(168).std()
    return X


def targets(c, h=HORIZON):
    c = pd.Series(c)
    lr = np.log(c).diff()
    fwd = np.log(c).shift(-h) - np.log(c)                    # forward return
    fvol = lr.shift(-h).rolling(h).std().shift(-(h - 1))     # forward realized vol
    return fwd, fvol


def shuffle(c, v, rng):
    """Destroy time-ordering, preserve the return distribution AND the return/volume
    pairing. This is the control arm: same statistics, no temporal structure."""
    r = np.diff(np.log(np.asarray(c, float)))
    idx = rng.permutation(len(r))
    return float(c[0]) * np.exp(np.concatenate([[0], np.cumsum(r[idx])])), \
           np.concatenate([[v[0]], np.asarray(v, float)[1:][idx]])


def evaluate(c, v, label):
    X = features(c, v)
    fwd, fvol = targets(c)
    y_dir = (fwd > 0).astype(int)
    y_vol = (fvol > fvol.median()).astype(int)

    ok = X.notna().all(axis=1) & fwd.notna() & fvol.notna()
    X, y_dir, y_vol = X[ok], y_dir[ok], y_vol[ok]
    n = len(X); cut = int(n * 0.70)                          # walk-forward, no shuffle
    if cut < 500 or n - cut < 200:
        return None

    out = {}
    for name, y in (("dir", y_dir), ("vol", y_vol)):
        m = lgb.LGBMClassifier(n_estimators=300, learning_rate=0.03, num_leaves=31,
                               min_child_samples=100, subsample=0.8, colsample_bytree=0.8,
                               random_state=SEED, verbose=-1)
        m.fit(X.iloc[:cut], y.iloc[:cut])
        p = m.predict_proba(X.iloc[cut:])[:, 1]
        yt = y.iloc[cut:]
        out[name] = roc_auc_score(yt, p) if yt.nunique() > 1 else np.nan
    return out


if __name__ == "__main__":
    rng = np.random.default_rng(SEED)
    data = {}
    for f in sorted(glob.glob(f"{FUT}/*-1h-futures.feather")):
        a = os.path.basename(f).split("_")[0]
        d = pd.read_feather(f).sort_values("date").drop_duplicates("date")
        if len(d) < 3000 or (d["volume"] * d["close"]).median() * 24 < 2e6:
            continue
        data[a] = d.reset_index(drop=True)

    print(f"{len(data)} liquid perps | 1h | predict {HORIZON}h ahead")
    print("out-of-sample = last 30% of each series, walk-forward (never shuffled in time)\n")
    print("=" * 88)
    print(f"{'':7s} | {'DIRECTION (which way)':^31s} | {'VOLATILITY (how much)':^31s}")
    print(f"{'pair':7s} | {'real AUC':>10s} {'shuffled':>10s} {'edge':>8s} |"
          f" {'real AUC':>10s} {'shuffled':>10s} {'edge':>8s}")
    print("-" * 88)

    R = {"dir": [], "vol": []}
    S = {"dir": [], "vol": []}
    for a, d in sorted(data.items()):
        c, v = d["close"].values.astype(float), d["volume"].values.astype(float)
        real = evaluate(c, v, a)
        cs, vs = shuffle(c, v, rng)
        ctrl = evaluate(cs, vs, a + "*")
        if real is None or ctrl is None:
            continue
        for k in ("dir", "vol"):
            R[k].append(real[k]); S[k].append(ctrl[k])
        print(f"{a:7s} | {real['dir']:10.4f} {ctrl['dir']:10.4f} {real['dir']-ctrl['dir']:+8.4f} |"
              f" {real['vol']:10.4f} {ctrl['vol']:10.4f} {real['vol']-ctrl['vol']:+8.4f}")

    print("-" * 88)
    print(f"{'MEAN':7s} | {np.mean(R['dir']):10.4f} {np.mean(S['dir']):10.4f}"
          f" {np.mean(R['dir'])-np.mean(S['dir']):+8.4f} |"
          f" {np.mean(R['vol']):10.4f} {np.mean(S['vol']):10.4f}"
          f" {np.mean(R['vol'])-np.mean(S['vol']):+8.4f}")
    print("\n  AUC 0.50 = coin flip.  AUC 1.00 = perfect.")
    print(f"  DIRECTION : real {np.mean(R['dir']):.4f} vs control {np.mean(S['dir']):.4f}")
    print(f"  VOLATILITY: real {np.mean(R['vol']):.4f} vs control {np.mean(S['vol']):.4f}")