Warrior_EA/research/learn_selection.py

"""Ask the model the answerable question: is the edge bigger than the cost RIGHT NOW?

Every neural net in this project so far has been asked to predict direction, and direction
is the one thing 23 years of data says is barely predictable. But the fade of retail
pin/inside-bar setups carries a MEASURED gross edge of ~0.139 R, and the cost of taking it
is not a constant - it is `spread / stop distance`, which varies bar to bar by a factor of
five across instruments, sessions and volatility regimes.

That makes trade SELECTION a well-posed supervised problem where direction never was:

    label    the realised R of the fade - what actually happened, in the unit that pays
    features everything knowable at the fill bar, including the cost itself
    goal     take the subset where edge > cost, skip the rest

The model is not being asked to find an edge. The edge is already measured. It is being
asked to spend it only where it survives - which is exactly what the per-cell table showed
a human doing by hand when it picked EURUSD over gold.

WHAT WOULD MAKE THIS SELF-DECEPTION, and is therefore controlled for:
  - Leakage. Every feature is computed from bars at or before the fill bar. The volume
    profile levels come from SESSIONS STRICTLY BEFORE the current one.
  - Shuffled splits. Forbidden here. Train / validate / test are chronological blocks, and
    the threshold is chosen on validation and applied unchanged to test.
  - Reporting the training fit. Only the test block is quoted.
  - The cost feature doing all the work. A model that learns nothing except "skip wide
    spreads" is still useful, but it is not a pattern - so the run reports what happens
    with the cost feature REMOVED, which separates the two.
"""
import numpy as np, sys, os, datetime as dt
sys.stdout.reconfigure(encoding='utf-8', errors='replace')
from test_retail import setups, triggered, race_px, load_bars, PIP
from vplevels import broker_day

SYMS = ('EURUSD', 'USDJPY', 'XAUUSD', 'SP500')
NAMES = ['cost', 'risk_atr', 'atr_rel', 'atr_pct', 'spread_pct', 'hour', 'dow',
         'is_pin', 'is_inside', 'is_engulf', 'dir', 'body', 'uwick', 'lwick',
         'ma_dist', 'ma_slope', 'don20', 'don50', 'ret5', 'ret20',
         'd_naked', 'd_vah', 'd_val', 'd_hvn', 'd_lvn', 'tickrate', 'rvol_rel', 'sym']


def _pct(x, n=500):
    """Rolling percentile rank of x within its own last n values. Causal."""
    out = np.full(len(x), 0.5)
    W = np.lib.stride_tricks.sliding_window_view
    if len(x) > n:
        w = W(x, n)[:-1]
        out[n:] = (w < x[n:, None]).mean(axis=1)
    return out


def vp_features(sym, t_ms):
    """Signed distance from close to each volume-profile level, in price units.

    Levels for broker-day d come only from days < d, so nothing here can see the session
    it is being used in.
    """
    try:
        from vplevels import load, day_index, session_levels, composite_nodes, naked_vpocs
        day, b, tk, dw, binsize = load(sym)
    except FileNotFoundError:
        return None
    days, cells = day_index(day, b, tk)
    sess = session_levels(days, cells)
    nodes = composite_nodes(days, cells, lookback=20)
    naked = naked_vpocs(sess)
    sday = {int(d): i for i, d in enumerate(sess[:, 0])}
    bd = broker_day(t_ms)
    return sess, nodes, naked, sday, binsize, bd


def build(sym, tf, H=200, path_tf='M5'):
    """One symbol/timeframe -> (X, y, times). y is the realised R of the FADE."""
    ev, o, h, l, c, spm, tick = setups(sym, tf)
    a1, I1 = load_bars(sym, tf)
    t1 = a1[:, I1['time']].astype(np.int64)
    a2, I2 = load_bars(sym, path_tf)
    ph, pl, pc = a2[:, I2['high']], a2[:, I2['low']], a2[:, I2['close']]
    pmap = np.searchsorted(a2[:, I2['time']], t1)
    HH = H * (12 if tf == 'H1' else 3)

    from test_cause_effect import atr_of
    atr = atr_of(h, l, c, 14); atr = np.concatenate([[atr[0]], atr[:-1]])
    ma = np.convolve(c, np.ones(20) / 20, mode='full')[:len(c)]; ma[:20] = c[:20]
    atrp = _pct(atr); spp = _pct(spm)
    rng_ = np.maximum(h - l, 1e-12)
    W = np.lib.stride_tricks.sliding_window_view
    d20 = np.full(len(c), 0.5); d50 = np.full(len(c), 0.5)
    for n, dst in ((20, d20), (50, d50)):
        if len(c) > n:
            hi = W(h, n).max(axis=1)[:-1]; lo = W(l, n).min(axis=1)[:-1]
            dst[n:] = (c[n:] - lo) / np.maximum(hi - lo, 1e-12)
    r5 = np.zeros(len(c)); r5[5:] = (c[5:] - c[:-5]) / np.maximum(atr[5:], 1e-12)
    r20 = np.zeros(len(c)); r20[20:] = (c[20:] - c[:-20]) / np.maximum(atr[20:], 1e-12)
    ticks = a1[:, I1['ticks']]; rvol = a1[:, I1['rvol']]

    vp = vp_features(sym, t1)
    X, Y, T = [], [], []
    for name, idx, d, ent, stp in ev:
        ok, fill = triggered(h, l, idx, d, ent)
        if ok.sum() < 100:
            continue
        i2, d2, e2, s2, sig = fill[ok], d[ok], ent[ok], stp[ok], idx[ok]
        risk0 = np.abs(e2 - s2)
        g = risk0 > 2 * spm[i2]
        i2, d2, e2, risk0, sig = (v[g] for v in (i2, d2, e2, risk0, sig))
        pi = np.clip(pmap[i2], 0, len(ph) - 1)
        keep = pi + HH < len(ph)
        i2, d2, e2, risk0, sig, pi = (v[keep] for v in (i2, d2, e2, risk0, sig, pi))
        if len(pi) < 200:
            continue
        sp = spm[i2]
        dd = -d2                                    # the fade
        r = race_px(ph, pl, pi, dd, e2 - dd * risk0, e2 + dd * risk0, HH)
        R = np.where(r > 0, 1.0, np.where(r < 0, -1.0, 0.0))
        un = r == 0
        if un.any():
            q = np.minimum(pi[un] + HH, len(pc) - 1)
            R[un] = (pc[q] - e2[un]) * dd[un] / risk0[un]
        R = R - sp / risk0

        #--- EVERY bar-derived feature must come from i2-1, not i2.
        #--- The order fills INTRABAR during bar i2, and the outcome race starts at the
        #--- first M5 bar inside i2 - so bar i2's close, tick count and realised volatility
        #--- are not knowable at entry. Reading them let the model see how the bar it
        #--- entered on turned out, and it scored +0.53 R on a held-out block. That is what
        #--- a lookahead looks like when it survives a clean chronological split: the split
        #--- was honest, the features were not.
        j2 = np.maximum(i2 - 1, 0)
        sp_known = spm[j2]                          # expected cost, not the realised one
        tt = t1[i2]
        hrs = np.array([dt.datetime.fromtimestamp(x / 1000, dt.UTC).hour for x in tt])
        dow = np.array([dt.datetime.fromtimestamp(x / 1000, dt.UTC).weekday() for x in tt])
        f = [sp_known / risk0, risk0 / atr[j2], atr[j2] / c[j2], atrp[j2], spp[j2],
             hrs, dow,
             np.full(len(i2), name == 'pin', float),
             np.full(len(i2), name == 'inside', float),
             np.full(len(i2), name == 'engulf', float),
             d2.astype(float),
             np.abs(c[sig] - o[sig]) / rng_[sig],
             (h[sig] - np.maximum(o[sig], c[sig])) / rng_[sig],
             (np.minimum(o[sig], c[sig]) - l[sig]) / rng_[sig],
             (c[j2] - ma[j2]) / np.maximum(atr[j2], 1e-12),
             (ma[j2] - ma[np.maximum(j2 - 5, 0)]) / np.maximum(atr[j2], 1e-12),
             d20[j2], d50[j2], r5[j2], r20[j2]]
        if vp:
            sess, nodes, naked, sday, binsize, bdall = vp
            bdi = bdall[i2]
            nk = np.zeros(len(i2)); vah = np.zeros(len(i2)); val = np.zeros(len(i2))
            hv = np.zeros(len(i2)); lv = np.zeros(len(i2))
            for q_ in range(len(i2)):
                dq = int(bdi[q_]); si = sday.get(dq)
                px = c[j2[q_]]; A = max(atr[j2[q_]], 1e-12)
                if si is not None and si >= 1:
                    vah[q_] = (sess[si - 1, 3] * binsize - px) / A
                    val[q_] = (sess[si - 1, 2] * binsize - px) / A
                nn = naked.get(dq)
                if nn is not None and len(nn):
                    nk[q_] = (nn[np.argmin(np.abs(nn * binsize - px))] * binsize - px) / A
                nd = nodes.get(dq)
                if nd is not None:
                    for arr, dst in ((nd[0], hv), (nd[1], lv)):
                        if len(arr):
                            dst[q_] = (arr[np.argmin(np.abs(arr * binsize - px))]
                                       * binsize - px) / A
            f += [nk, vah, val, hv, lv]
        else:
            f += [np.zeros(len(i2))] * 5
        f += [ticks[j2] / np.maximum(np.median(ticks), 1),
              rvol[j2] / np.maximum(atr[j2] ** 2, 1e-12),
              np.full(len(i2), SYMS.index(sym), float)]
        X.append(np.column_stack(f)); Y.append(R); T.append(tt)
    if not X:
        return None
    return np.vstack(X), np.concatenate(Y), np.concatenate(T)


def run(tfs=('M15', 'H1'), syms=SYMS, drop_cost=False, seed=0):
    Xs, Ys, Ts = [], [], []
    for tf in tfs:
        for s in syms:
            r = build(s, tf)
            if r:
                Xs.append(r[0]); Ys.append(r[1]); Ts.append(r[2])
                print(f"  {s} {tf}: {len(r[1]):,} fade trades, mean R {r[1].mean():+.4f}")
    X = np.vstack(Xs); y = np.concatenate(Ys); t = np.concatenate(Ts)
    o = np.argsort(t, kind='stable')               # chronological, pooled across symbols
    X, y, t = X[o], y[o], t[o]
    if drop_cost:
        X = X[:, 1:]
    n = len(y)
    a, b = int(0.60 * n), int(0.80 * n)
    U = lambda ms: dt.datetime.fromtimestamp(ms / 1000, dt.UTC).date()
    print(f"\n  {n:,} trades   train {U(t[0])}..{U(t[a-1])}   "
          f"val ..{U(t[b-1])}   TEST {U(t[b])}..{U(t[-1])}")
    print(f"  baseline mean R: train {y[:a].mean():+.4f}  val {y[a:b].mean():+.4f}  "
          f"TEST {y[b:].mean():+.4f}")

    from sklearn.ensemble import HistGradientBoostingRegressor
    m = HistGradientBoostingRegressor(max_iter=400, learning_rate=0.05, max_depth=6,
                                      min_samples_leaf=200, l2_regularization=1.0,
                                      random_state=seed)
    m.fit(X[:a], y[:a])
    pv, pt = m.predict(X[a:b]), m.predict(X[b:])
    #--- threshold chosen on VALIDATION only, then frozen
    best = None
    for q in np.arange(0.0, 0.95, 0.05):
        thr = np.quantile(pv, q)
        sel = pv >= thr
        if sel.sum() < 200:
            continue
        sc = y[a:b][sel].mean()
        if best is None or sc > best[1]:
            best = (thr, sc, q)
    thr, vsc, q = best
    sel = pt >= thr
    kept = y[b:][sel]
    se = kept.std(ddof=1) / np.sqrt(max(len(kept), 1))
    print(f"\n  threshold from validation: keep top {100*(1-q):.0f}%  (val mean R {vsc:+.4f})")
    print(f"  TEST: kept {len(kept):,}/{len(pt):,} ({100*sel.mean():.0f}%)   "
          f"mean R {kept.mean():+.4f} +/- {se:.4f}   t {kept.mean()/max(se,1e-9):+.2f}")
    print(f"  TEST all trades for comparison:              mean R {y[b:].mean():+.4f}")

    #--- The t above is not believable as stated: with an 8-day horizon these trades overlap
    #--- heavily, so thousands of them share the same price path and the effective sample is
    #--- far smaller than the count. Two honest checks.
    tt, yy = t[b:][sel], y[b:][sel]
    yrs = np.array([dt.datetime.fromtimestamp(x / 1000, dt.UTC).year for x in tt])
    print("  by year:  " + "  ".join(
        f"{u}:{yy[yrs == u].mean():+.3f}(n={int((yrs == u).sum())})"
        for u in np.unique(yrs)))
    pos = sum(1 for u in np.unique(yrs) if yy[yrs == u].mean() > 0)
    print(f"  years positive: {pos}/{len(np.unique(yrs))}")

    #--- non-overlapping subset: keep a trade only if the previous kept one has expired
    gap = 200 * 3600 * 1000            # ~the H1 horizon, in ms; conservative for M15 too
    keep, busy = [], -1
    for q_ in range(len(tt)):
        if tt[q_] < busy:
            continue
        keep.append(q_); busy = tt[q_] + gap
    keep = np.array(keep, int)
    z = yy[keep]
    se2 = z.std(ddof=1) / np.sqrt(max(len(z), 1))
    print(f"  NON-OVERLAPPING: {len(z)} independent trades   mean R {z.mean():+.4f}"
          f" +/- {se2:.4f}   t {z.mean()/max(se2,1e-9):+.2f}")
    return m, X, y, t, b, sel


if __name__ == '__main__':
    print("=== SELECTING FADE TRADES: can a model spend the edge only where it survives? ===\n")
    run()
    print("\n=== SAME, WITH THE COST FEATURE REMOVED ===")
    print("    if this collapses, the model learned 'avoid wide spreads' and nothing else -")
    print("    useful, but not a pattern.\n")
    run(drop_cost=True)