backtest: real CuPy GPU backend (kernels are bit-identical CPU twins); books byte-identical on both backends (2026-09-18)

CuPy from pip had never compiled a kernel on the dev box (no CUDA_PATH for the
splayed CUDA-12 wheel layout); the probe swallowed the error and every book to
date was built on the NumPy path under a "GPU backtester" name.

- backtest/gpu_kernels.py (new): lazily probed, never-silent BACKEND state
  (CUDA_PATH pointed at the nvidia-cuda-runtime-cu12 header wheel; the probe
  compiles a kernel and keeps the failure reason), --backend auto|gpu|cpu and
  FB2_BACKEND, and three CUDA kernels compiled with --fmad=false: ewm_batch
  (the _gpu_ewm_alpha recurrence, one thread per series), roll_mean_var
  (pandas 3.0.2 roll_mean/roll_var replica, Kahan + Welford + unstable
  recompute), resolve_trades (trailing and static exits, one thread per trade).
- gpu_backtest.py: precompute_indicators and simulate_trades dispatch to the
  GPU twins or to the untouched NumPy paths; caches are host NumPy and device
  memory is released before the replay workers fork; the Backend: line names
  the device or the fallback reason.
- test_all.py Module 36 (7 tests): a visible CUDA device must never fall back
  silently; kernels vs their CPU references on synthetic edge cases; backend
  selection; import guards (cupy banned from the live roots). Suite 536.
- CLAUDE.md, instructions_for_backtest.txt, research note with the evidence.

Verification: all five reference MD5s reproduced on the GPU and with
--backend cpu (ten runs, GPU cache fallbacks 0); precompute GPU == CPU on all
42 cached windows (2.34 M bars); suite green with GPU, with
CUDA_VISIBLE_DEVICES="" and with FB2_BACKEND=cpu; --workers 1 and --record
pairs equal. The per-bar replay stays the live Python signal code on CPU
workers by design; wall time is roughly unchanged (precompute ~3x faster).
This commit is contained in:
Harwinder Sekhon
2026-09-18 10:39:55 -04:00
parent ad14ea4fed
commit 02589cbb88
6 changed files with 1487 additions and 70 deletions
+307 -60
View File
@@ -19,7 +19,9 @@ Key properties:
- Trade simulation is full candle walk on the primary timeframe using
High/Low; if both SL and TP1 hit in the same bar, SL is assumed first
(conservative).
- Ships with CuPy fallback to NumPy so the script runs anywhere.
- Ships with CuPy fallback to NumPy so the script runs anywhere. The CUDA
kernels (backtest/gpu_kernels.py) are bit-for-bit replicas of the NumPy
path, so both backends produce the same book; --backend cpu forces NumPy.
Usage:
python backtest/gpu_backtest.py --pair XAUUSD --days 365
@@ -64,36 +66,22 @@ from engine import IndicatorSignal, TechnicalAnalysis, TechnicalEngine # noqa:
from signal_aggregator import SignalAggregator, SignalReport, TradeSetup # noqa: E402
import position_sizing # noqa: E402
from backtest import replay_record # noqa: E402
from backtest import gpu_kernels as _gk # noqa: E402
# ─────────────────────────────────────────────────────────────────────
# GPU shim — CuPy if available AND CUDA runtime works, else NumPy.
# An importable CuPy can still fail at first kernel launch if CUDA libs
# are missing (containers, headless boxes), so verify with a tiny op.
# GPU backend — backtest/gpu_kernels.py owns the CuPy state (probed
# lazily and never silently: a fallback carries its reason into the
# "Backend:" line) and the CUDA kernels. Every kernel is a bit-for-bit
# replica of the NumPy/pandas code kept in this file as the CPU path, so
# `--backend cpu` and the GPU produce the same book (proof: the five
# reference MD5s, CLAUDE.md section 9). The forked replay workers never
# touch the device: every precomputed series is moved to host first.
# ─────────────────────────────────────────────────────────────────────
def _detect_gpu():
try:
import cupy as cp_mod # type: ignore
# Force runtime init: a kernel launch will surface missing libs.
_ = cp_mod.asarray(np.zeros(2, dtype=np.float64))
_ = cp_mod.cumsum(_).get()
return cp_mod, True
except Exception:
return np, False
BACKEND = _gk.BACKEND
cp, _GPU_AVAILABLE = _detect_gpu()
def _to_gpu(arr: np.ndarray):
"""Move a numpy array to GPU (or no-op if CuPy unavailable)."""
return cp.asarray(arr) if _GPU_AVAILABLE else arr
def _to_cpu(arr) -> np.ndarray:
"""Pull array back to host."""
if _GPU_AVAILABLE and hasattr(arr, "get"):
return arr.get()
return np.asarray(arr)
def _gpu_on() -> bool:
return BACKEND.available
logger = logging.getLogger("gpu_backtest")
@@ -215,18 +203,29 @@ class IndicatorCache:
def precompute_indicators(
df: pd.DataFrame, tech_cfg, gpu_log: bool = False
) -> IndicatorCache:
"""Compute every causal indicator series once on the full df."""
"""Compute every causal indicator series once on the full df.
Dispatches to the CUDA kernels when the backend is available and to
the NumPy/pandas code otherwise. The two are bit-identical (Module 36
and the 42-window check in backtest/research/gpu_backend_2026-09-18.md).
"""
if gpu_log:
logger.info(
f"Pre-computing indicators on {len(df)} bars "
f"({'CuPy GPU' if _gpu_on() else 'NumPy CPU'} backend)"
)
if _gpu_on():
return _precompute_indicators_gpu(df, tech_cfg)
return _precompute_indicators_cpu(df, tech_cfg)
def _precompute_indicators_cpu(df: pd.DataFrame, tech_cfg) -> IndicatorCache:
"""The NumPy/pandas path — the reference the CUDA twin replicates."""
close = df["Close"].to_numpy(dtype=np.float64)
high = df["High"].to_numpy(dtype=np.float64)
low = df["Low"].to_numpy(dtype=np.float64)
n = len(close)
if gpu_log:
logger.info(
f"Pre-computing indicators on {n} bars "
f"({'GPU' if _GPU_AVAILABLE else 'CPU'} backend)"
)
# ── EMAs (use ewm span semantics, identical to pandas .ewm(span=, adjust=False).mean()) ─
ema_fast = _gpu_ewm_span(close, tech_cfg.ema_fast)
ema_medium = _gpu_ewm_span(close, tech_cfg.ema_medium)
@@ -388,6 +387,187 @@ def precompute_indicators(
)
def _assert_host_cache(cache: IndicatorCache) -> None:
"""Every array the forked workers will read must be plain NumPy."""
for name, val in vars(cache).items():
if name in ("df", "index_pos", "rsi_series", "macd_hist_series") or val is None:
continue
if type(val) is not np.ndarray:
raise TypeError(f"IndicatorCache.{name} is {type(val).__name__}, not np.ndarray")
def _precompute_indicators_gpu(df: pd.DataFrame, tech_cfg) -> IndicatorCache:
"""CUDA twin of `_precompute_indicators_cpu`.
Same expressions in the same order (single-op CuPy ufuncs cannot fuse
into an FMA and double add/mul/div/sqrt are correctly rounded on both
sides); the recurrences and the pandas rolling statistics run in the
kernels of backtest/gpu_kernels.py, one thread per series. Everything
comes back to the host once at the end: the replay workers read NumPy.
"""
cp = BACKEND.cp
nan = float("nan")
close_h = df["Close"].to_numpy(dtype=np.float64)
n = len(close_h)
close = cp.asarray(close_h)
high = cp.asarray(df["High"].to_numpy(dtype=np.float64))
low = cp.asarray(df["Low"].to_numpy(dtype=np.float64))
# ── RSI inputs (delta[0] is NaN; NaN > 0 is False so gain[0] = 0.0) ──
delta = cp.empty(n, dtype=cp.float64)
delta[0] = nan
delta[1:] = close[1:] - close[:-1]
gain = cp.where(delta > 0, delta, 0.0)
loss = cp.where(delta < 0, -delta, 0.0)
# ── True range (fmax ignores NaN, like the CPU path) ──
prev_close = cp.empty(n, dtype=cp.float64)
prev_close[0] = nan
prev_close[1:] = close[:-1]
tr1 = high - low
tr2 = cp.abs(high - prev_close)
tr3 = cp.abs(low - prev_close)
tr2[0] = nan
tr3[0] = nan
tr = cp.fmax(cp.fmax(tr1, tr2), tr3)
# ── Directional movement ──
plus_dm_raw = cp.empty(n, dtype=cp.float64)
minus_dm_raw = cp.empty(n, dtype=cp.float64)
plus_dm_raw[0] = nan
minus_dm_raw[0] = nan
plus_dm_raw[1:] = high[1:] - high[:-1]
minus_dm_raw[1:] = -(low[1:] - low[:-1])
cond_plus = (plus_dm_raw > minus_dm_raw) & (plus_dm_raw > 0)
plus_dm = cp.where(cond_plus, plus_dm_raw, 0.0)
cond_minus = (minus_dm_raw > plus_dm_raw) & (minus_dm_raw > 0)
minus_dm = cp.where(cond_minus, minus_dm_raw, 0.0)
# ── Stage 1: every recurrence that depends only on price, one launch ──
def span_alpha(span: int) -> float:
return 2.0 / (span + 1.0) # identical to _gpu_ewm_span
rsi_p = tech_cfg.rsi_period
atr_p = tech_cfg.atr_period
adx_p = tech_cfg.adx_period
kc_period = getattr(tech_cfg, 'kc_period', 20)
kc_atr_period = getattr(tech_cfg, 'kc_atr_period', 10)
kc_multiplier = getattr(tech_cfg, 'kc_multiplier', 1.5)
have_trend = n >= tech_cfg.ema_trend
jobs = [
(close, span_alpha(tech_cfg.ema_fast), 0),
(close, span_alpha(tech_cfg.ema_medium), 0),
(close, span_alpha(tech_cfg.ema_slow), 0),
(close, span_alpha(tech_cfg.macd_fast), 0),
(close, span_alpha(tech_cfg.macd_slow), 0),
(close, span_alpha(kc_period), 0),
(gain, 1.0 / rsi_p, rsi_p),
(loss, 1.0 / rsi_p, rsi_p),
(tr, 1.0 / atr_p, atr_p),
(tr, 1.0 / adx_p, adx_p),
(plus_dm, 1.0 / adx_p, adx_p),
(minus_dm, 1.0 / adx_p, adx_p),
(tr, 1.0 / kc_atr_period, kc_atr_period),
]
if have_trend:
jobs.append((close, span_alpha(tech_cfg.ema_trend), 0))
out = _gk.ewm_batch_dev([j[0] for j in jobs], [j[1] for j in jobs],
[j[2] for j in jobs])
(ema_fast, ema_medium, ema_slow, ema_macd_fast, ema_macd_slow, kc_mid,
avg_gain, avg_loss, atr, atr_adx, plus_dm_smooth, minus_dm_smooth,
kc_atr) = out[:13]
ema_trend = out[13] if have_trend else None
# ── RSI ──
avg_loss_safe = cp.where(avg_loss == 0, 1e-10, avg_loss)
rs = avg_gain / avg_loss_safe
rsi = 100.0 - (100.0 / (1.0 + rs))
# ── MACD line; DI / DX ──
macd_line = ema_macd_fast - ema_macd_slow
atr_adx_safe = cp.where(atr_adx == 0, nan, atr_adx)
plus_di = 100.0 * plus_dm_smooth / atr_adx_safe
minus_di = 100.0 * minus_dm_smooth / atr_adx_safe
di_sum = plus_di + minus_di
di_sum_safe = cp.where(di_sum == 0, nan, di_sum)
dx = 100.0 * cp.abs(plus_di - minus_di) / di_sum_safe
# ── Stage 2: recurrences on stage-1 outputs, one launch ──
signal_line, adx = _gk.ewm_batch_dev(
[macd_line, dx], [span_alpha(tech_cfg.macd_signal), 1.0 / adx_p], [0, adx_p])
macd_hist = macd_line - signal_line
# ── Bollinger: the pandas rolling replica (mean + var, then zsqrt) ──
(bb_mid,), (bb_var,) = _gk.roll_mean_var_dev([close], tech_cfg.bb_period, ddof=0)
bb_std = _gk.roll_std_from_var_dev(bb_var)
bb_upper = bb_mid + tech_cfg.bb_std * bb_std
bb_lower = bb_mid - tech_cfg.bb_std * bb_std
bb_range = bb_upper - bb_lower
bb_range_safe = cp.where(bb_range == 0, nan, bb_range)
bb_position = (close - bb_lower) / bb_range_safe
bb_position = cp.where(cp.isnan(bb_position), 0.5, bb_position)
bb_mid_safe = cp.where(bb_mid == 0, nan, bb_mid)
bb_band_width = (bb_upper - bb_lower) / bb_mid_safe
# 20-bar mean of the band width: the CPU path uses the literal 20.
(bb_avg_width,), _unused = _gk.roll_mean_var_dev([bb_band_width], 20, ddof=0)
# ── ADX slope ──
adx_slope_lookback = getattr(tech_cfg, 'adx_slope_lookback', 5)
adx_slope = cp.full(n, nan, dtype=cp.float64)
if n > adx_slope_lookback:
adx_slope[adx_slope_lookback:] = adx[adx_slope_lookback:] - adx[:-adx_slope_lookback]
# ── ROC ──
roc_period = getattr(tech_cfg, 'roc_period', 10)
roc = cp.full(n, nan, dtype=cp.float64)
if n > roc_period:
prev_close = close[:-roc_period]
prev_safe = cp.where(prev_close == 0, nan, prev_close)
roc[roc_period:] = ((close[roc_period:] - prev_close) / prev_safe) * 100.0
# ── Keltner ──
kc_upper = kc_mid + kc_multiplier * kc_atr
kc_lower = kc_mid - kc_multiplier * kc_atr
kc_squeeze = ((bb_lower > kc_lower) & (bb_upper < kc_upper)).astype(cp.float64)
# ── One transfer to the host; the workers read NumPy only ──
def h(a):
return a.get()
rsi_h = h(rsi)
macd_hist_h = h(macd_hist)
cache = IndicatorCache(
df=df,
index_pos={ts: i for i, ts in enumerate(df.index)},
ema_fast=h(ema_fast),
ema_medium=h(ema_medium),
ema_slow=h(ema_slow),
ema_trend=h(ema_trend) if ema_trend is not None else None,
rsi=rsi_h,
rsi_series=pd.Series(rsi_h, index=df.index, name="rsi"),
macd_hist=macd_hist_h,
macd_hist_series=pd.Series(macd_hist_h, index=df.index, name="macd_hist"),
bb_mid=h(bb_mid),
bb_upper=h(bb_upper),
bb_lower=h(bb_lower),
bb_position=h(bb_position),
bb_band_width=h(bb_band_width),
bb_avg_width=h(bb_avg_width),
atr=h(atr),
adx=h(adx),
plus_di=h(plus_di),
minus_di=h(minus_di),
adx_slope=h(adx_slope),
roc=h(roc),
kc_mid=h(kc_mid),
kc_upper=h(kc_upper),
kc_lower=h(kc_lower),
kc_squeeze=h(kc_squeeze),
)
_assert_host_cache(cache)
return cache
# ═════════════════════════════════════════════════════════════════════
# GPU-AWARE TECHNICAL ENGINE (subclasses live engine, lookup-only)
# ═════════════════════════════════════════════════════════════════════
@@ -1330,6 +1510,15 @@ def simulate_trades(
closes_arr[i, :n] = sim_close[entry_pos:end_pos]
exit_times.append(list(sim_index[entry_pos:end_pos]))
if _gpu_on():
# One CUDA thread per trade, both exit modes (bit-identical twin of
# the two CPU paths below; Module 36 and the reference MD5s).
_resolve_trades_gpu(
trades, opens_arr, highs_arr, lows_arr, closes_arr,
valid_lens, exit_times, spread, trail_arm_r, trail_dist_r,
)
return
if trail_arm_r > 0:
_simulate_trailing(
trades, opens_arr, highs_arr, lows_arr, closes_arr,
@@ -1337,42 +1526,33 @@ def simulate_trades(
)
return
# GPU vectorized: for each trade, scan its window for first SL/TP hit.
# CPU static path, vectorized: for each trade, scan its window for the
# first SL/TP hit.
sls = np.array([t.sl for t in trades], dtype=np.float64)[:, None]
tps = np.array([t.tp1 for t in trades], dtype=np.float64)[:, None]
is_buy = np.array([t.direction == "BUY" for t in trades], dtype=bool)[:, None]
h_g = _to_gpu(highs_arr)
l_g = _to_gpu(lows_arr)
c_g = _to_gpu(closes_arr)
sls_g = _to_gpu(sls)
tps_g = _to_gpu(tps)
is_buy_g = _to_gpu(is_buy)
# Hit masks. For BUY: SL when low<=sl, TP when high>=tp.
# For SELL: SL when high>=sl, TP when low<=tp.
sl_hit_buy = l_g <= sls_g
tp_hit_buy = h_g >= tps_g
sl_hit_sell = h_g >= sls_g
tp_hit_sell = l_g <= tps_g
sl_hit_buy = lows_arr <= sls
tp_hit_buy = highs_arr >= tps
sl_hit_sell = highs_arr >= sls
tp_hit_sell = lows_arr <= tps
sl_hit = cp.where(is_buy_g, sl_hit_buy, sl_hit_sell)
tp_hit = cp.where(is_buy_g, tp_hit_buy, tp_hit_sell)
sl_hit = np.where(is_buy, sl_hit_buy, sl_hit_sell)
tp_hit = np.where(is_buy, tp_hit_buy, tp_hit_sell)
# Mask out NaN bars (no data past valid_lens).
valid_mask = ~cp.isnan(h_g)
valid_mask = ~np.isnan(highs_arr)
sl_hit = sl_hit & valid_mask
tp_hit = tp_hit & valid_mask
# First-hit indices (or max_bars if never).
idx_template = cp.arange(max_bars)[None, :]
sl_idx = cp.where(sl_hit, idx_template, max_bars)
tp_idx = cp.where(tp_hit, idx_template, max_bars)
first_sl = cp.min(sl_idx, axis=1)
first_tp = cp.min(tp_idx, axis=1)
sl_first_arr = _to_cpu(first_sl)
tp_first_arr = _to_cpu(first_tp)
idx_template = np.arange(max_bars)[None, :]
sl_idx = np.where(sl_hit, idx_template, max_bars)
tp_idx = np.where(tp_hit, idx_template, max_bars)
sl_first_arr = np.min(sl_idx, axis=1)
tp_first_arr = np.min(tp_idx, axis=1)
# Resolve each trade and compute realized MFE/MAE *up to exit bar*
# (a trade closed at bar 5 cannot have meaningful excursions at bar 12).
@@ -1447,6 +1627,60 @@ def simulate_trades(
t.mae_r = float(max(0.0, mae))
def _resolve_trades_gpu(
trades: List[Trade],
opens_arr: np.ndarray,
highs_arr: np.ndarray,
lows_arr: np.ndarray,
closes_arr: np.ndarray,
valid_lens: np.ndarray,
exit_times: List[List[pd.Timestamp]],
spread: float,
trail_arm_r: float,
trail_dist_r: float,
) -> None:
"""CUDA twin of `_simulate_trailing` (trail_arm_r > 0) and of the static
SL/TP1 branch of `simulate_trades`: the kernel walks every trade's window
in its own thread and returns exit index/price/reason, bars held, R and
the excursions; this host loop only writes them onto the Trade objects,
with the same 'not set' cases as the CPU code."""
trailing = trail_arm_r > 0
entries = np.array([t.entry_price for t in trades], dtype=np.float64)
sls = np.array([t.sl for t in trades], dtype=np.float64)
tps = np.array([t.tp1 for t in trades], dtype=np.float64)
is_buy = np.array([t.direction == "BUY" for t in trades], dtype=np.int8)
arms = np.array([t.trail_arm_r if t.trail_arm_r > 0 else trail_arm_r
for t in trades], dtype=np.float64)
dists = np.array([t.trail_dist_r if t.trail_dist_r > 0 else trail_dist_r
for t in trades], dtype=np.float64)
res = _gk.resolve_trades(
opens_arr, highs_arr, lows_arr, closes_arr, valid_lens,
entries, sls, tps, is_buy, arms, dists, spread, trailing)
for i, t in enumerate(trades):
code = int(res["reason"][i])
if code == _gk.REASON_NO_DATA:
t.exit_reason = "NO_DATA"
continue
bars = int(res["bars_held"][i])
exit_idx = int(res["exit_idx"][i])
if code == _gk.REASON_NO_DATA_AFTER_TIMEOUT:
# Static path: the TIMEOUT bookkeeping ran before the NaN close
# was seen, so bars_held and exit_time are kept and the price
# is None (never a NaN R in the book).
t.exit_reason = "NO_DATA"
t.bars_held = bars
t.exit_price = None
t.exit_time = exit_times[i][exit_idx]
continue
t.exit_reason = _gk.REASON_NAMES[code]
t.bars_held = bars
t.exit_price = float(res["exit_price"][i])
t.exit_time = exit_times[i][exit_idx]
t.r_multiple = float(res["r"][i])
t.mfe_r = float(res["mfe"][i])
t.mae_r = float(res["mae"][i])
# ═════════════════════════════════════════════════════════════════════
# BACKTEST RUNNER
# ═════════════════════════════════════════════════════════════════════
@@ -1878,7 +2112,12 @@ def run_backtest(
tf.value: precompute_indicators(dfs_full[tf], platform_cfg.technical, gpu_log=True)
for tf in timeframes
}
logger.info(f"Indicator pre-compute took {time.perf_counter() - t0:.2f}s")
logger.info(
f"Indicator pre-compute took {time.perf_counter() - t0:.2f}s "
f"({'CuPy GPU' if _gpu_on() else 'NumPy CPU'})")
# The replay workers are forked below and must never touch the device:
# every cache array is already host NumPy; hand pooled memory back too.
BACKEND.release_device_memory()
# ── Inject GPU engine into a fresh aggregator ──
aggregator = SignalAggregator(platform_cfg)
@@ -1998,7 +2237,7 @@ def run_backtest(
trail_dist_r=platform_cfg.signals.trail_dist_r)
logger.info(
f"Trade simulation took {time.perf_counter() - t_sim:.2f}s "
f"(GPU={'yes' if _GPU_AVAILABLE else 'no'})"
f"({'CuPy GPU' if _gpu_on() else 'NumPy CPU'})"
)
if bt_cfg.record_dir:
@@ -2009,7 +2248,7 @@ def run_backtest(
n_primary_bars=len(primary_df), n_bars_in_window=n_bars,
skipped=skipped_reasons, n_candidates=len(candidates),
workers=workers,
backend="CuPy GPU" if _GPU_AVAILABLE else "NumPy CPU fallback",
backend=BACKEND.describe(),
cache_files=cache_names, replay_seconds=replay_seconds,
cache_fallbacks=cache_fallbacks),
sim_df=sim_df, bar_chunks=bar_chunks or None)
@@ -2421,6 +2660,12 @@ def _parse_args(argv: Optional[List[str]] = None) -> argparse.Namespace:
help="Print full per-trade reasoning (indicators, signals, agreement, why entered)")
p.add_argument("--workers", type=int, default=0,
help="Replay worker processes (0=auto: cpu-2 capped at 12, 1=sequential)")
p.add_argument("--backend", choices=("auto", "gpu", "cpu"), default=None,
help="Compute backend for the indicator pre-compute and the "
"trade resolution: auto = CuPy when CUDA works, else "
"NumPy (the reason is logged); gpu = fail loudly without "
"CUDA; cpu = force the NumPy twins (parity runs). The "
"env var FB2_BACKEND sets the default.")
p.add_argument("--reasoning-samples", type=int, default=3,
help="If --detail not set, print this many sample trades' full reasoning")
p.add_argument("--log-level", default="INFO")
@@ -2486,6 +2731,8 @@ def main(argv: Optional[List[str]] = None) -> int:
level=args.log_level.upper(),
format="%(asctime)s [%(levelname)s] %(name)s: %(message)s",
)
if args.backend:
BACKEND.set_mode(args.backend)
pair = _resolve_pair(args.pair)
primary_tf = _resolve_tf(args.primary_tf)
@@ -2592,7 +2839,7 @@ def main(argv: Optional[List[str]] = None) -> int:
)
logger.info(
f"Backend: {'CuPy GPU' if _GPU_AVAILABLE else 'NumPy CPU fallback'} | "
f"Backend: {BACKEND.describe()} | "
f"Pair: {PAIR_INFO[pair]['display']} | Window: {start.date()} → {end.date()} | "
f"Confirm TFs: {','.join(tf.value for tf in confirm_tfs)} | Spread: {spread}"
)