Warrior_EA/research/ticks_to_bars.py

"""Stream SQX tick files into bars carrying tick-derived microstructure features.

One pass, bounded memory. A full symbol is ~458M ticks (~15 GB as arrays), so nothing can
hold the raw stream; this reduces 23 years of EURUSD ticks to ~2.4M M5 bars that every
downstream test can load instantly.

WHAT IS AND IS NOT COMPUTABLE HERE
----------------------------------
The feed carries (time, bid, ask, volume). It does NOT carry trade direction: SQX's own
record has one volume field, and MT5's TICK_FLAG_BUY/SELL are empty on FX anyway. So true
signed order flow does not exist in this data and is not synthesised here under a
flattering name.

What IS available is quote dynamics, which is the standard substitute when trade labels
are absent:

  tick rule     - sign of each mid-price change; up-ticks vs down-ticks over the bar.
                  The classic Lee-Ready fallback classifier.
  OFI counts    - SIGNED top-of-book revisions: bid_up, bid_dn, ask_up, ask_dn. A bid
                  ticking up and an ask ticking down both mean buy-side pressure, and a
                  bare "the bid changed" counter cannot tell them apart. This is
                  order-flow imbalance in the Cont/Kukanov/Stoikov sense, in its
                  EVENT-COUNT form - the feed has no sizes, so it cannot be size-weighted.
                  Temper expectations accordingly: OFI is well established as a
                  CONTEMPORANEOUS explainer of price change, and its predictive power
                  decays within seconds. At M5 with multi-hour horizons the prior should
                  be that it explains the bar it is measured in and not the next one.
  arrival rate  - ticks per bar, and the dispersion of inter-tick gaps. Urgency, not size.
  realised vol  - sum of squared mid returns within the bar; a far better volatility
                  estimate than the bar range, and only obtainable from ticks.
  spread stats  - mean and max within the bar, in price units.

All of these are UNSIGNED except the tick rule and quote asymmetry, which are the only two
columns here that can point a direction.

Bars are stamped by the OPEN of their interval and computed only from ticks inside it, so
nothing in a bar's features depends on a tick after it closes.
"""
import numpy as np, sys, os, time, datetime as dt
sys.stdout.reconfigure(encoding='utf-8', errors='replace')
from sqx import decode_iter, calibrate_decimals, NO_VOLUME, VOLUME_CONSTANT

COLUMNS = ['time', 'open', 'high', 'low', 'close', 'ticks', 'upticks', 'downticks',
           'bid_up', 'bid_dn', 'ask_up', 'ask_dn', 'spread_mean', 'spread_max',
           'rvol', 'volume', 'gap_mean', 'gap_max']
#--- column groups used when merging a bar that straddles a batch boundary
_I = {c: k for k, c in enumerate(COLUMNS)}
SUM_IDX = tuple(_I[c] for c in ('ticks', 'upticks', 'downticks', 'bid_up', 'bid_dn',
                                'ask_up', 'ask_dn', 'spread_mean', 'rvol', 'volume',
                                'gap_mean'))
MAX_IDX = tuple(_I[c] for c in ('spread_max', 'gap_max'))
MEAN_IDX = tuple(_I[c] for c in ('spread_mean', 'gap_mean'))   # accumulated as sums, divided at the end


def merge(a, b):
    """Combine two partials describing the SAME bucket. a precedes b.
    Index-driven off COLUMNS rather than positional, so adding a feature cannot
    silently mis-merge one."""
    out = list(a)
    out[_I['high']] = max(a[_I['high']], b[_I['high']])
    out[_I['low']] = min(a[_I['low']], b[_I['low']])
    out[_I['close']] = b[_I['close']]
    for i in SUM_IDX:
        out[i] = a[i] + b[i]
    for i in MAX_IDX:
        out[i] = max(a[i], b[i])
    return tuple(out)


def finalise(arr):
    """Sort, merge duplicate-timestamp rows, then turn accumulated sums into means.

    Called ONCE on the fully assembled array. Chunk workers must not divide their own
    means: a bar straddling a chunk boundary arrives as two partials, and mean-of-means
    weighted wrong is exactly the kind of small silent error that survives every sanity
    check.

    The sort is needed because the raw feed is not perfectly time-ordered - EURUSD carries
    2 backward steps in 513M ticks, both under an hour, both in the 2003-2006 era. Bucketing
    is by ABSOLUTE timestamp, so every tick still lands in its correct bar no matter what
    order it arrives in; the only symptom is that a bucket can be emitted twice, out of
    order. A stable sort followed by the duplicate merge below puts those back together
    exactly. The one residual imprecision is that `open`/`close` for such a bar come from
    decode order rather than true time order - at most a couple of bars in 1.7M."""
    if not len(arr):
        return arr
    t = arr[:, _I['time']]
    if (np.diff(t) < 0).any():
        arr = arr[np.argsort(t, kind='stable')]
        t = arr[:, _I['time']]
    if (np.diff(t) == 0).any():
        starts = np.concatenate(([0], np.flatnonzero(np.diff(t)) + 1))
        ends = np.concatenate((starts[1:] - 1, [len(t) - 1]))
        out = arr[starts].copy()
        out[:, _I['high']] = np.maximum.reduceat(arr[:, _I['high']], starts)
        out[:, _I['low']] = np.minimum.reduceat(arr[:, _I['low']], starts)
        out[:, _I['close']] = arr[ends, _I['close']]
        for i in SUM_IDX:
            out[:, i] = np.add.reduceat(arr[:, i], starts)
        for i in MAX_IDX:
            out[:, i] = np.maximum.reduceat(arr[:, i], starts)
        arr = out
    n = np.maximum(arr[:, _I['ticks']], 1.0)
    for i in MEAN_IDX:
        arr[:, i] /= n
    return arr


def _chunk(args):
    """One byte range of one file -> partial bars with means still UNDIVIDED."""
    path, scale, step, start, stop, volume_constant, max_records = args
    return _scan(path, scale, step, volume_constant, max_records, start, stop)[0]


def _scan(path, scale, step, volume_constant, max_records=None, start=None, stop=None,
          progress_every=0, t0=None):
    """Decode a byte range and aggregate it to bars. Returns (array, tick_count)."""
    rows = []
    pend = None                       # carried partial bar as a plain tuple
    prev_mid = prev_bid = prev_ask = None
    prev_t = None
    total = 0
    t0 = t0 or time.time()

    for ts, ai, bi, vi in decode_iter(path, max_records=max_records, start=start, stop=stop):
        ask = ai / scale
        bid = bi / scale
        mid = (ask + bid) * 0.5
        vol = np.where(vi == NO_VOLUME, 0.0, vi / volume_constant)
        spread = ask - bid
        # per-tick deltas, seeded from the previous batch's last tick so the first tick of
        # a batch is not silently treated as having no predecessor
        d = np.empty(len(mid))
        d[0] = 0.0 if prev_mid is None else mid[0] - prev_mid
        d[1:] = np.diff(mid)
        # SIGNED quote changes, not just "did it move". Order-flow imbalance in the
        # Cont/Kukanov/Stoikov sense is built from the DIRECTION of top-of-book revisions:
        # a bid ticking UP and an ask ticking DOWN both mean buy-side pressure, and the two
        # are indistinguishable from a bare "the bid changed" counter. Sizes are absent
        # from this feed, so this is the event-count form of OFI rather than the
        # size-weighted one - see the module docstring for what that costs.
        db = np.empty(len(bid))
        da = np.empty(len(ask))
        db[0] = 0.0 if prev_bid is None else bid[0] - prev_bid
        da[0] = 0.0 if prev_ask is None else ask[0] - prev_ask
        db[1:] = np.diff(bid)
        da[1:] = np.diff(ask)
        bid_up = (db > 0).astype(np.float64)
        bid_dn = (db < 0).astype(np.float64)
        ask_up = (da > 0).astype(np.float64)
        ask_dn = (da < 0).astype(np.float64)
        gap = np.empty(len(ts))
        gap[0] = 0.0 if prev_t is None else (ts[0] - prev_t) / 1000.0
        gap[1:] = np.diff(ts) / 1000.0
        np.maximum(gap, 0.0, out=gap)

        bucket = (ts // step) * step
        starts = np.concatenate(([0], np.flatnonzero(np.diff(bucket)) + 1))
        ends = np.concatenate((starts[1:] - 1, [len(ts) - 1]))
        cnt = (ends - starts + 1).astype(np.float64)

        part = np.column_stack([
            bucket[starts].astype(np.float64),        # time
            mid[starts],                              # open
            np.maximum.reduceat(mid, starts),         # high
            np.minimum.reduceat(mid, starts),         # low
            mid[ends],                                # close
            cnt,                                      # ticks
            np.add.reduceat((d > 0).astype(np.float64), starts),   # upticks
            np.add.reduceat((d < 0).astype(np.float64), starts),   # downticks
            np.add.reduceat(bid_up, starts),          # bid ticked up
            np.add.reduceat(bid_dn, starts),          # bid ticked down
            np.add.reduceat(ask_up, starts),          # ask ticked up
            np.add.reduceat(ask_dn, starts),          # ask ticked down
            np.add.reduceat(spread, starts),          # spread SUM (divided at the end)
            np.maximum.reduceat(spread, starts),      # spread max
            np.add.reduceat(d * d, starts),           # realised variance
            np.add.reduceat(vol, starts),             # volume
            np.add.reduceat(gap, starts),             # gap SUM
            np.maximum.reduceat(gap, starts),         # gap max
        ])
        lst = [tuple(r) for r in part]
        if pend is not None:
            if lst and lst[0][0] == pend[0]:
                lst[0] = merge(pend, lst[0])
            else:
                rows.append(pend)
        pend = lst.pop() if lst else pend
        rows.extend(lst)

        prev_mid = mid[-1]; prev_bid = bid[-1]; prev_ask = ask[-1]; prev_t = ts[-1]
        total += len(ts)
        if progress_every and total % progress_every < len(ts):
            el = time.time() - t0
            print(f"  {total/1e6:.0f}M ticks  {len(rows):,} bars  "
                  f"{dt.datetime.fromtimestamp(ts[-1]/1000, dt.UTC):%Y-%m-%d}  "
                  f"{total/el/1e6:.2f}M/s", flush=True)
    if pend is not None:
        rows.append(pend)
    return np.array(rows, dtype=np.float64), total


def build(path, ref, tf_seconds=300, decimals=None, max_records=None,
          volume_constant=VOLUME_CONSTANT, progress_every=20000000, workers=None):
    """ref = (time_s, close) reference series used only to calibrate the price scale.

    The per-record decode loop is pure Python and cannot be vectorised - record LENGTH
    depends on the config nibbles, so you cannot know where record k+1 starts without
    parsing record k. It therefore saturates exactly one core, which is why this ran at
    ~10% CPU on a 12-core box. But the format is randomly seekable: every BLOCK_LENGTH
    records SQX restates all four fields as absolute int64s, so byte ranges starting at
    block headers decode independently with no shared history. That is the whole trick.

    The only cost is that the first tick of each chunk has no predecessor, so its per-tick
    deltas (tick direction, quote revisions, inter-tick gap) count as zero. That is
    workers-1 ticks out of ~458M, i.e. ~2e-8 of the sample - stated here rather than left
    for someone to find.
    """
    from sqx import decode, decode_iter, split_offsets, find_block
    # --- calibrate on a sample first; the scale must be right before a long pass
    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:
            # The head of the tick file can predate the MT5 reference series entirely -
            # XAUUSD ticks start 2003-05-05 but its H1 export starts 2004-06-11, so the
            # first 200k ticks overlap it by nothing. Recent data always overlaps, so fall
            # back to a sample near the END of the file. Failing here rather than guessing
            # is deliberate: a wrong scale crashes nothing and silently rescales every
            # ATR-normalised feature downstream.
            size = os.path.getsize(path)
            print("  head sample predates the reference series - resampling")
            #--- 2M ticks, because the sample must span >=50 reference HOURS, not >=50
            #--- ticks: 200k ticks of modern gold is about nine hours and calibration
            #--- rightly refused it. Midpoint first (always inside the reference range),
            #--- then the tail.
            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
    step = tf_seconds * 1000
    t0 = time.time()

    if workers is None:
        workers = max(1, (os.cpu_count() or 4) - 1)
    if max_records or workers <= 1:
        arr, total = _scan(path, scale, step, volume_constant, max_records,
                           progress_every=progress_every, t0=t0)
    else:
        ranges = split_offsets(path, workers)
        print(f"  {len(ranges)} chunks x ~{os.path.getsize(path)/len(ranges)/1e9:.2f} GB "
              f"across {workers} workers", flush=True)
        import multiprocessing as mp
        args = [(path, scale, step, s, e, volume_constant, None) for s, e in ranges]
        with mp.Pool(len(ranges)) as pool:
            parts = []
            for k, part in enumerate(pool.imap(_chunk, args)):
                parts.append(part)
                el = time.time() - t0
                print(f"  chunk {k+1}/{len(ranges)} done  {len(part):,} bars  "
                      f"{el/60:.1f} min elapsed", flush=True)
        parts = [p for p in parts if len(p)]
        arr = np.concatenate(parts) if parts else np.empty((0, len(COLUMNS)))
        total = int(arr[:, _I['ticks']].sum()) if len(arr) else 0
        #--- imap preserves submission order and ranges are in file order, so the
        #--- concatenation should be time-sorted. It is not quite, because the FEED is not
        #--- perfectly ordered. Distinguish the two causes by magnitude rather than by
        #--- presence: a genuine split bug misplaces whole chunks and puts a large fraction
        #--- of rows out of order, while feed noise puts a handful. finalise() sorts either
        #--- way, but only one of them is acceptable to continue past.
        if len(arr):
            ooo = int((np.diff(arr[:, _I['time']]) < 0).sum())
            if ooo > max(16, len(arr) // 1000):
                raise RuntimeError(f"{ooo:,}/{len(arr):,} bars out of time order - "
                                   "that is a chunking bug, not feed noise")
            if ooo:
                print(f"  {ooo} out-of-order bar(s) from feed noise - sorting")

    before = len(arr)
    arr = finalise(arr)
    if before != len(arr):
        print(f"  merged {before - len(arr)} bar(s) straddling chunk boundaries")
    print(f"  DONE {total:,} ticks -> {len(arr):,} bars in {(time.time()-t0)/60:.1f} min")
    return arr


if __name__ == '__main__':
    from kit import load_rates
    D = 'c:/Users/admin/Documents/Workspaces/Market Data/'
    OUT = 'c:/Users/admin/Documents/Workspaces/Market Data/bars/'
    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')]
    tf = int(sys.argv[1]) if len(sys.argv) > 1 else 300
    only = sys.argv[2] if len(sys.argv) > 2 else None
    for fn, sym in JOBS:
        if only and sym != only:
            continue
        out = f"{OUT}{sym}_M{tf//60}_ticks.npz"
        if os.path.exists(out):
            print(f"{sym}: {out} exists, skipping")
            continue
        print(f"\n=== {sym} ({os.path.getsize(D+fn)/1e9:.1f} GB) -> M{tf//60} ===", flush=True)
        rt, ro, rh, rl, rc, rv, rs = load_rates(sym, 16385)
        arr = build(D + fn, (rt, rc), tf_seconds=tf)
        np.savez_compressed(out, bars=arr, columns=np.array(COLUMNS))
        print(f"  saved {out}")