"""Is the EA's single volume feature leaving information on the table?

The EA feeds exactly one volume value per bar:

    volumeChangeRatio = (v[i] - v[i-1]) / v[i-1]      (clamped +/-5)

That is the first difference. It cannot express the level relative to a baseline, it is
blind to the intraday volume profile - the strongest and most reliable pattern in FX tick
volume - and it cannot distinguish "high volume, big range" (continuation) from "high
volume, small range" (absorption), which mean opposite things.

This measures whether richer representations carry more information ABOUT THE LABEL than
the shipped one. Deliberately not a "fit a model and look for edge" test: this session has
twice produced a 4-sigma result that was a lookahead, so the question here is the narrower,
harder-to-fool one - does the feature share mutual information with the outcome, beyond
what its own finite-sample bias and the label's autocorrelation would produce anyway.

Null: BLOCK permutation of the labels, block length = the barrier horizon. Barrier labels
on adjacent bars share almost all of their outcome window, so a free shuffle destroys that
autocorrelation and produces a null that is far too tight - every feature then looks
significant. Blocks preserve it.

Bias: MI of independent variables is not 0 in finite samples but approximately
(bins-1)(classes-1)/(2n). With 8 bins and 3 classes that is 7/n. Reported, not subtracted -
the permutation null already absorbs it, and quoting both is what makes the number readable.
"""
import numpy as np, sys, time
sys.stdout.reconfigure(encoding='utf-8', errors='replace')
from kit import load_rates, atr, sma, barrier_vec

BINS = 8
NPERM = 500


def rank_bin(x, bins=BINS):
    """Equal-frequency (rank) binning, so the estimate does not depend on the feature's
    marginal shape - only on how the label sorts across it. Done ONCE per feature: the
    bin assignment does not change when the LABELS are permuted, and re-deriving it inside
    the permutation loop was what made the first version of this unrunnable."""
    n = len(x)
    r = np.argsort(np.argsort(x))
    return (r * bins // n).clip(0, bins - 1).astype(np.int64)


def mi_binned(xb, y, bins=BINS):
    """MI from a pre-binned feature. Contingency table via one bincount - O(n), no masks."""
    n = len(y)
    if n < 100:
        return 0.0
    joint = np.bincount(xb * 3 + y, minlength=bins * 3).reshape(bins, 3).astype(float) / n
    px = joint.sum(axis=1, keepdims=True)
    py = joint.sum(axis=0, keepdims=True)
    nz = joint > 0
    return float((joint[nz] * np.log(joint[nz] / (px @ py)[nz])).sum())


def block_perm_null_multi(xbs, y, horizon, nperm=NPERM, seed=0):
    """Shuffle BLOCKS of labels, preserving within-block order.

    One permutation is applied to EVERY feature before drawing the next, so all features
    see the identical null draws. That makes their nulls comparable (and is what a
    family-wise max-statistic would need), at the cost of correlating them - which is fine
    here because each feature is judged against its own null, not against the others.
    """
    n = len(y)
    blk = max(int(horizon), 1)
    nb = (n + blk - 1) // blk
    rng = np.random.default_rng(seed)
    out = np.empty((len(xbs), nperm))
    starts = np.arange(nb) * blk
    for b in range(nperm):
        order = rng.permutation(nb)
        # gather block start offsets in shuffled order, then expand to indices
        perm = np.concatenate([np.arange(starts[o], min(starts[o] + blk, n)) for o in order])[:n]
        yp = y[perm]
        for f in range(len(xbs)):
            out[f, b] = mi_binned(xbs[f][:len(yp)], yp)
    return out


def session_z(v, hours, lookback_days=20):
    """Volume z-score against the SAME hour-of-day's own recent history.

    Causal by construction: for each bar, statistics come only from previous occurrences
    of that hour. This is the feature the EA has no way to express - the network currently
    has to infer the intraday profile from the hour sin/cos features and combine it with a
    raw volume ratio itself.
    """
    n = len(v)
    out = np.zeros(n)
    for h in range(24):
        m = np.nonzero(hours == h)[0]
        if len(m) < 5:
            continue
        vh = v[m].astype(float)
        # trailing mean/sd over the last `lookback_days` occurrences of this hour, shifted
        # by one so the current bar never contributes to its own baseline
        cs = np.concatenate([[0.0], np.cumsum(vh)])
        cs2 = np.concatenate([[0.0], np.cumsum(vh * vh)])
        for k in range(len(m)):
            lo = max(0, k - lookback_days)
            cnt = k - lo
            if cnt < 5:
                continue
            s = cs[k] - cs[lo]
            s2 = cs2[k] - cs2[lo]
            mu = s / cnt
            var = max(s2 / cnt - mu * mu, 1e-12)
            out[m[k]] = (vh[k] - mu) / np.sqrt(var)
    return np.clip(out, -5, 5)


def build_features(o, h, l, c, v, t):
    a = atr(h, l, c, 14)
    a = np.where(a > 0, a, np.nan)
    hours = np.array([(int(x) // 3600) % 24 for x in t])
    prev = np.roll(v, 1).astype(float); prev[0] = v[0]
    base = sma(v.astype(float), 50)
    rng_atr = (h - l) / a
    vlevel = np.where(base > 0, v / base, 1.0)

    F, names = [], []
    F.append(np.clip(np.where(prev > 0, (v - prev) / prev, 0.0), -5, 5)); names.append("volChange (SHIPPED)")
    F.append(np.clip(vlevel, 0, 5));                                     names.append("volLevel  v/sma50")
    F.append(session_z(v, hours));                                       names.append("volSessionZ")
    # absorption: range delivered per unit of volume. Low = lots of activity, little travel.
    F.append(np.clip(np.where(vlevel > 0.05, rng_atr / vlevel, 0.0), 0, 5)); names.append("absorption rng/vol")
    # the interaction the shipped feature cannot express at all
    F.append(np.clip(vlevel * rng_atr, 0, 5));                           names.append("vol x range")
    return [np.nan_to_num(f) for f in F], names


def run(sym, tf, sl_m, tp_m, H):
    t, o, h, l, c, v, spr = load_rates(sym, tf)
    a = atr(h, l, c, 14)
    tick = np.nanmin(np.abs(np.diff(np.unique(np.round(c, 8)))))
    sp = np.nanmedian(spr) * tick
    if not np.isfinite(sp):
        sp = 0.0
    lab, valid = barrier_vec(h, l, c, np.where(a > 0, a, np.nan), sl_m, tp_m, H, sp)
    F, names = build_features(o, h, l, c, v, t)
    m = valid & np.isfinite(a) & (a > 0)
    m[:200] = False
    m[-(H + 2):] = False
    y = lab[m]
    n = int(m.sum())
    print(f"\n=== {sym}  SL{sl_m}:TP{tp_m} H={H}  n={n}  "
          f"labels Buy {100*(y==0).mean():.1f}% Sell {100*(y==1).mean():.1f}% Neu {100*(y==2).mean():.1f}% ===")
    print(f"  finite-sample MI bias ~ 7/n = {7.0/n:.6f} nats")
    xbs = [rank_bin(f[m]) for f in F]
    obs = [mi_binned(xb, y) for xb in xbs]
    nulls = block_perm_null_multi(xbs, y, H)
    print(f"{'feature':<22}{'MI (nats)':>12}{'null mean':>12}{'null p95':>11}{'excess':>10}{'p':>8}")
    for k, nm in enumerate(names):
        null = nulls[k]
        p = (1 + int((null >= obs[k]).sum())) / (NPERM + 1)   # Phipson & Smyth
        star = ' *' if p < 0.05 else ''
        print(f"{nm:<22}{obs[k]:>12.6f}{null.mean():>12.6f}{np.quantile(null,0.95):>11.6f}"
              f"{obs[k]-null.mean():>+10.6f}{p:>8.3f}{star}", flush=True)


if __name__ == '__main__':
    t0 = time.time()
    for sym in ('EURUSD', 'USDJPY', 'XAUUSD'):
        for (s, p, H) in [(2, 3, 96), (1, 2, 64)]:
            try:
                run(sym, 16385, s, p, H)
            except Exception as ex:
                print(f"{sym} {s}:{p} FAILED {ex}")
    print(f"\ntotal {time.time()-t0:.0f}s")