Warrior_EA/research/profiles.py

"""Build true volume-at-price profiles from the raw tick stream.

Everything this project has tested so far reduced the market to a BAR: OHLC plus some
aggregate. A volume profile is the one description that a bar cannot carry, because it is
volume distributed ACROSS PRICE inside the period rather than summed over it. Villahermosa's
"Wyckoff 2.0" is built entirely on that distribution - VPOC, value area, high/low volume
nodes - so testing any of its claims requires computing it properly first.

WHAT IS BEING COUNTED, HONESTLY
-------------------------------
The feed has no trade prints and no trade direction (see ticks_to_bars.py). So this is a
count of QUOTE UPDATES at each price level, not contracts traded. That distinction matters
and is not papered over:

  - It is closest in spirit to Steidlmayer's TPO / Market Profile, which is time-at-price
    and needs no volume at all. Several of the claims under test (the 80% rule, value-area
    acceptance) are ORIGINALLY TPO claims, so this is the right measurement for them, not
    a compromise.
  - For the volume-specific claims (VPOC as institutional inventory) it is a proxy. Quote
    activity and traded volume correlate strongly but are not identical.
  - A price level where the market sat quietly for an hour accumulates fewer counts than
    one where it thrashed for ten minutes. That is a real difference from a time profile
    and it is the reason both weightings are computed here: `ticks` and `dwell`
    (milliseconds spent at the level), so any result can be checked against the other.

BINNING
-------
Bins are on a FIXED ABSOLUTE GRID (bin = round(price / size)), never relative to the day's
range. Two profiles can only be compared - and a node can only be "the same level" across
sessions - if their bins line up. A per-day adaptive grid would silently destroy exactly
the cross-session persistence being tested.

DAY BOUNDARY
------------
Broker day, i.e. UTC+2 (measured, see the stop-run notes). The daily bar the5ers shows
rolls at 22:00 UTC and every "previous session value area" claim is about THAT session,
not about a calendar day in UTC.
"""
import numpy as np, sys, os, time, datetime as dt
sys.stdout.reconfigure(encoding='utf-8', errors='replace')
from sqx import decode_iter, split_offsets, calibrate_decimals, decode, find_block

BARS = 'c:/Users/admin/Documents/Workspaces/Market Data/bars/'
OUT = 'c:/Users/admin/Documents/Workspaces/Market Data/profiles/'
BROKER_OFFSET_MS = 2 * 3600 * 1000        # broker = UTC + 2
DAY_MS = 86400 * 1000

#--- Bin size per instrument, chosen so a typical DAY spans O(1000) bins: fine enough that
#--- a VPOC is a price rather than a zone, coarse enough that counts are not shot noise.
#--- Explicit rather than derived, for the same reason the decimals table is explicit.
BIN = {'EURUSD': 0.00001,    # 0.1 pip   - day range ~0.008 ->  ~800 bins
       'USDJPY': 0.001,      # 0.1 pip   - day range ~0.9   ->  ~900 bins
       'XAUUSD': 0.01,       # 1 cent    - day range ~20    -> ~2000 bins
       'SP500':  0.10}       # 0.1 index - day range ~50    ->  ~500 bins


def _scan(args):
    """One byte range -> (key, ticks, dwell) where key packs (day, bin).

    Returned already reduced by np.unique, so a chunk covering 40M ticks comes back as a
    few hundred thousand rows instead of 40M. The reduction is exact - counts are summed,
    never sampled.
    """
    path, scale, binsize, start, stop = args
    keys, tks, dws = [], [], []
    for ts, ai, bi, _vi in decode_iter(path, start=start, stop=stop):
        mid = (ai + bi) * (0.5 / scale)
        day = (ts + BROKER_OFFSET_MS) // DAY_MS
        b = np.rint(mid / binsize).astype(np.int64)
        #--- pack: bin can be ~300k for gold, day ~20k. 2^32 is ample headroom for both,
        #--- and an int64 key keeps the unique() fast. Bins are non-negative for every
        #--- instrument here (prices are positive), asserted below.
        key = day * (1 << 32) + b
        #--- dwell = time this tick's price stood before the next quote arrived. The last
        #--- tick of a chunk has no successor inside the chunk; it gets the chunk's median
        #--- gap rather than zero, which is worth ~1 tick in 40M and avoids a systematic
        #--- (if tiny) bias against the final price of every chunk.
        gap = np.empty(len(ts), np.float64)
        gap[:-1] = np.diff(ts)
        gap[-1] = np.median(gap[:-1]) if len(ts) > 1 else 0.0
        np.clip(gap, 0, 60000, out=gap)   # a weekend gap is not dwell at a price
        if (b < 0).any():
            raise ValueError("negative price bin - check the decimals calibration")
        k, inv = np.unique(key, return_inverse=True)
        keys.append(k)
        tks.append(np.bincount(inv, minlength=len(k)).astype(np.float64))
        dws.append(np.bincount(inv, weights=gap, minlength=len(k)))
    if not keys:
        return (np.empty(0, np.int64), np.empty(0), np.empty(0))
    return _reduce(np.concatenate(keys), np.concatenate(tks), np.concatenate(dws))


def _reduce(key, tk, dw):
    k, inv = np.unique(key, return_inverse=True)
    return k, np.bincount(inv, weights=tk, minlength=len(k)), \
           np.bincount(inv, weights=dw, minlength=len(k))


def build(path, sym, ref, decimals=None, workers=None):
    binsize = BIN[sym]
    if decimals is None:
        sts, sai, _, _ = decode(path, max_records=200000)
        try:
            decimals, err = calibrate_decimals(sts, sai, ref[0], ref[1], verbose=False)
        except ValueError:
            size = os.path.getsize(path)
            decimals = None
            for frac in (0.5, 0.75, 0.9):
                with open(path, 'rb') as fh:
                    off = find_block(fh, int(size * frac), size)
                if off is None:
                    continue
                it = decode_iter(path, batch=2000000, start=off)
                sts, sai, _, _ = next(it); it.close()
                try:
                    decimals, err = calibrate_decimals(sts, sai, ref[0], ref[1],
                                                       verbose=False)
                    break
                except ValueError:
                    continue
            if decimals is None:
                raise
        print(f"  scale: decimals={decimals} (median rel.err {err:.4%})")
    scale = 10.0 ** decimals
    workers = workers or max(1, (os.cpu_count() or 4) - 1)
    ranges = split_offsets(path, workers)
    t0 = time.time()
    print(f"  {len(ranges)} chunks across {workers} workers, bin={binsize}", flush=True)
    import multiprocessing as mp
    args = [(path, scale, binsize, s, e) for s, e in ranges]
    parts = []
    with mp.Pool(len(ranges)) as pool:
        for i, p in enumerate(pool.imap(_scan, args)):
            parts.append(p)
            print(f"  chunk {i+1}/{len(ranges)}  {len(p[0]):,} cells  "
                  f"{(time.time()-t0)/60:.1f} min", flush=True)
    key, tk, dw = _reduce(np.concatenate([p[0] for p in parts]),
                          np.concatenate([p[1] for p in parts]),
                          np.concatenate([p[2] for p in parts]))
    day = (key >> 32).astype(np.int64)
    b = (key & ((1 << 32) - 1)).astype(np.int64)
    print(f"  {int(tk.sum()):,} ticks -> {len(key):,} (day,price) cells over "
          f"{len(np.unique(day)):,} days in {(time.time()-t0)/60:.1f} min")
    return day, b, tk, dw, binsize


#----------------------------------------------------------------------------------------
# Profile statistics
#----------------------------------------------------------------------------------------
def value_area(bins, w, frac=0.70):
    """Classic Market Profile value area: start at the POC and repeatedly annex whichever
    NEIGHBOURING PAIR of levels holds more weight, until `frac` of the total is enclosed.

    Not a percentile and not mean +/- sigma. Those give a different (usually wider) region
    on any skewed profile, and skew is the whole point of the P/b shapes. Villahermosa
    quotes 68.2% by analogy with one standard deviation; the Market Profile convention that
    every platform implements is 70%. The difference is immaterial to the tests here and
    70% is used so the numbers are comparable with anything else the user reads.
    """
    n = len(bins)
    if n == 0:
        return None
    poc = int(np.argmax(w))
    #--- prefix sums: the walk asks for a running total O(n) times, and recomputing the
    #--- enclosed sum each step turns a 350-bin day into 60k adds for nothing.
    cs = np.concatenate(([0.0], np.cumsum(w)))
    seg = lambda a, b: cs[b + 1] - cs[a]      # inclusive [a, b]
    tot = cs[-1]
    need = frac * tot
    lo = hi = poc
    got = w[poc]
    while got < need and (lo > 0 or hi < n - 1):
        #--- pairs, per the standard algorithm, so the area does not creep one side at a time
        dn = seg(max(lo - 2, 0), lo - 1) if lo > 0 else -1.0
        up = seg(hi + 1, min(hi + 2, n - 1)) if hi < n - 1 else -1.0
        if up >= dn:
            if hi >= n - 1:
                break
            hi = min(hi + 2, n - 1)
        else:
            if lo <= 0:
                break
            lo = max(lo - 2, 0)
        got = seg(lo, hi)
    return bins[poc], bins[lo], bins[hi]


def per_day(day, b, w):
    """-> (days, vpoc, val, vah, lo, hi, total) with bins as INTEGER bin indices.

    Gaps matter: a price level with zero activity inside the day's range must be present as
    a zero, otherwise the value-area walk annexes across a hole as if it were adjacent, and
    low-volume nodes - the entire basis of the stop-placement claim - become invisible.
    """
    order = np.lexsort((b, day))
    day, b, w = day[order], b[order], w[order]
    starts = np.concatenate(([0], np.flatnonzero(np.diff(day)) + 1))
    ends = np.concatenate((starts[1:], [len(day)]))
    D, P, L, H, LO, HI, T = [], [], [], [], [], [], []
    for s, e in zip(starts, ends):
        bb, ww = b[s:e], w[s:e]
        full = np.arange(bb[0], bb[-1] + 1)
        dense = np.zeros(len(full))
        dense[bb - bb[0]] = ww
        va = value_area(full, dense)
        D.append(day[s]); P.append(va[0]); L.append(va[1]); H.append(va[2])
        LO.append(bb[0]); HI.append(bb[-1]); T.append(ww.sum())
    return (np.array(D), np.array(P), np.array(L), np.array(H),
            np.array(LO), np.array(HI), np.array(T))


if __name__ == '__main__':
    from kit import load_rates
    D = 'c:/Users/admin/Documents/Workspaces/Market Data/'
    os.makedirs(OUT, exist_ok=True)
    JOBS = [('EURUSD_tick_the5ers_TICK.dat', 'EURUSD'),
            ('USDJPY_tick_the5ers_TICK.dat', 'USDJPY'),
            ('XAUUSD_tick_the5ers_TICK.dat', 'XAUUSD'),
            ('USA500IDXUSD_tick_the5ers_TICK.dat', 'SP500')]
    only = sys.argv[1] if len(sys.argv) > 1 else None
    for fn, sym in JOBS:
        if only and sym != only:
            continue
        out = f"{OUT}{sym}_vp.npz"
        if os.path.exists(out):
            print(f"{sym}: exists, skipping"); continue
        print(f"\n=== {sym} volume profile ({os.path.getsize(D+fn)/1e9:.1f} GB) ===",
              flush=True)
        rt, ro, rh, rl, rc, rv, rs = load_rates(sym, 16385)
        day, b, tk, dw, binsize = build(D + fn, sym, (rt, rc))
        np.savez_compressed(out, day=day, bin=b, ticks=tk, dwell=dw,
                            binsize=np.array([binsize]))
        print(f"  saved {out}")