Warrior_EA/research/test_mi_control.py

"""Is the EA's MI positive control failing because the ESTIMATOR is broken, or because its
own premise is wrong?

The EA prints, on every topology:

    "MI positive control - label 24 bars away ... scores 0.00307 nats against the ~0.00189
     noise floor. WARNING - the estimator FAILED to detect an association that must be
     there. Every mutual-information figure above is void."

That warning voids the barrier-geometry scan, which is the instrument that decides WHICH
target is worth training on. So it has to be resolved before another training cycle, and
there are two very different explanations:

  A. the estimator (or the row pairing) is broken, and the real MI is large
  B. two triple-barrier labels 24-48 bars apart genuinely share very little, the control's
     premise "must be strongly associated" is false, and the warning is a false alarm that
     is voiding perfectly good measurements

Code reading narrowed it but cannot decide it. `BuildMiSample` shifts BOTH ends of its bar
range by |offset| and recomputes the stride from that range, then the control pairs the two
builds BY INDEX - so the pair is actually 48 bars apart, not 24, and the stride can differ
between the two builds (it did not in this run: 19 for both). Neither of those turns a large
association into 0.003 nats on its own.

So: reproduce the label on validated data and measure the answer directly.

The estimator is reimplemented EXACTLY as `FeatureColumnMI` does it - rank into 8 equal-count
bins, plug-in MI on the 8x3 contingency table - and reported beside an exact discrete MI on
the 3x3 table. If the two agree, the estimator is sound and the answer is B.
"""
import numpy as np, sys
sys.stdout.reconfigure(encoding='utf-8', errors='replace')
import book, fills

MI_BINS = 8


def triple_barrier(h, l, c, atr, sl_mult, tp_mult, horizon):
    """0=Buy, 1=Sell, 2=Neutral. Long wins if +tp is touched before -sl, and vice versa.

    Matches the EA's convention: a bar containing both barriers counts as the STOP, which is
    why the two directions are resolved independently and a tie falls through to Neutral.
    """
    n = len(c)
    lab = np.full(n, 2, np.int8)
    ok = np.zeros(n, bool)
    INF = np.iinfo(np.int32).max
    lTp = c + tp_mult * atr; lSl = c - sl_mult * atr
    sTp = c - tp_mult * atr; sSl = c + sl_mult * atr
    CH = max(1, 4_000_000 // max(horizon, 1))
    for s in range(0, n, CH):
        e = min(s + CH, n - horizon)
        if e <= s:
            break
        idx = np.arange(s, e)
        w = idx[:, None] + np.arange(1, horizon + 1)[None, :]
        wh = h[w]; wl = l[w]

        def first(m):
            a = m.any(axis=1)
            return np.where(a, m.argmax(axis=1), INF)
        ltp = first(wh >= lTp[idx, None]); lsl = first(wl <= lSl[idx, None])
        stp = first(wl <= sTp[idx, None]); ssl = first(wh >= sSl[idx, None])
        lw = ltp < lsl; sw = stp < ssl
        seg = np.full(len(idx), 2, np.int8)
        seg[lw & ~sw] = 0
        seg[sw & ~lw] = 1
        lab[s:e] = seg
        ok[s:e] = np.isfinite(atr[idx]) & (atr[idx] > 0)
    return lab, ok


def mi_exact(a, b):
    """Exact plug-in MI of two discrete variables, in nats."""
    ka, kb = int(a.max()) + 1, int(b.max()) + 1
    j = np.bincount(a.astype(int) * kb + b.astype(int), minlength=ka * kb).reshape(ka, kb)
    n = j.sum()
    if n == 0:
        return 0.0
    p = j / n
    pa = p.sum(1, keepdims=True); pb = p.sum(0, keepdims=True)
    with np.errstate(divide='ignore', invalid='ignore'):
        t = p * np.log(p / (pa * pb))
    return float(np.nansum(np.where(p > 0, t, 0.0)))


def mi_ea(vals, labels):
    """The EA's estimator, reproduced: rank into MI_BINS equal-count bins, then plug-in MI.

    `rank` counts strictly-smaller elements, so ties share a bin - which for a 3-valued
    column means the three values land in three distinct bins, and no information is lost.
    """
    n = len(vals)
    order = np.argsort(vals, kind='stable')
    sv = vals[order]
    rank = np.searchsorted(sv, vals, side='left')
    bx = np.minimum((rank * MI_BINS // n), MI_BINS - 1)
    return mi_exact(bx.astype(int), labels.astype(int))


if __name__ == '__main__':
    sym, tf = 'SP500', 'H1'
    bk = fills.Book(sym)
    f = book.frame(sym, tf, bk)
    atr = f.atr(14)
    #--- the shipped geometry, read off the EA's own barrier-geometry scan line: 2:4, h=96
    lab, ok = triple_barrier(f.h, f.l, f.c, atr, 2.0, 4.0, 96)
    v = ok & (lab >= 0)
    sh = np.bincount(lab[v], minlength=3) / v.sum()
    print(f"=== {sym} {tf} triple-barrier 2:4 h=96 ===")
    print(f"  reproduced label shares  Buy {sh[0]:.1%}  Sell {sh[1]:.1%}  Neutral {sh[2]:.1%}")
    print(f"  EA reported              Buy 33.4%  Sell 29.7%  Neutral 36.9%")
    H = -np.sum(sh * np.log(np.maximum(sh, 1e-12)))
    print(f"  label entropy {H:.3f} nats (EA: 1.094)\n")

    #--- the EA samples every `stride` bars, so replicate that too: the control is computed on
    #--- a strided subsample, not on every bar, and a strided sample of an autocorrelated
    #--- series is exactly where an association can quietly disappear
    print("  MI between the label and the label k bars away, in nats")
    print(f"  {'k':>5}{'n':>9}{'exact 3x3':>12}{'EA 8-bin':>11}{'% of H':>9}   note")
    idx = np.nonzero(v)[0]
    for k, note in ((1, 'adjacent bar'), (24, 'the control CLAIMS this'),
                    (48, 'what it actually measures'), (96, 'one full horizon'),
                    (192, 'two horizons - must be ~0')):
        a = idx[idx + k < len(lab)]
        b = a + k
        m = v[b]
        a, b = a[m], b[m]
        if len(a) < 500:
            continue
        e = mi_exact(lab[a], lab[b])
        g = mi_ea(lab[b].astype(float), lab[a])
        print(f"  {k:>5}{len(a):>9,}{e:>12.5f}{g:>11.5f}{100*e/H:>8.1f}%   {note}")

    print("\n  same, on the EA's STRIDED subsample (stride 19, ~2000 rows):")
    print(f"  {'k':>5}{'n':>9}{'exact 3x3':>12}{'EA 8-bin':>11}{'% of H':>9}")
    sub = idx[::19][:2009]
    for k in (24, 48):
        a = sub[sub + k < len(lab)]
        b = a + k
        m = v[b]
        a, b = a[m], b[m]
        e = mi_exact(lab[a], lab[b])
        g = mi_ea(lab[b].astype(float), lab[a])
        print(f"  {k:>5}{len(a):>9,}{e:>12.5f}{g:>11.5f}{100*e/H:>8.1f}%")
    print("\n  If these land near 0.003 nats the control's premise is false and the WARNING is")
    print("  a false alarm voiding good measurements. If they land far above, the pairing or")
    print("  the estimator in the EA is broken.")