ForexBot2/backtest/research/sizing_model.py

"""Stage-D stake-sizing model driver (2026-09 campaign).

Walk-forward fit + evaluation of a per-signal stake multiplier on a trade
book exported by ``gpu_backtest.py --export-features`` (63 base columns +
``stake_multiplier``/``r_sized`` + 65 ``f_*`` decision-time features).

The protocol implemented here is the one pre-registered in
``backtest/research/sizing_campaign_2026-09.md`` and is NOT to be varied to
rescue a result:

* expanding walk-forward, 12-month burn-in, calendar-quarter refits, a
  ``purge_bars × bar_minutes`` gap between the end of train and the start of
  the test quarter, standardisation fit on train only;
* models: ridge on winsorised R (primary), L2 logistic on ``R > 0``
  (secondary), a one-feature ``f_is_buy`` simplicity baseline, an
  exhaustion-block-only ablation, and an optional lightgbm upper bound;
* size map ``m = round_0.05(clip(normaliser × (0.5 + clip((z−lo)/(hi−lo),0,1)),
  0.5, 1.5))`` with ``lo``/``hi``/``normaliser`` calibrated on TRAIN
  predictions only;
* a label-shuffle null (training labels permuted, test labels untouched);
* six acceptance clauses, printed PASS/FAIL individually.

Usage:
    python backtest/research/sizing_model.py --book BOOK.csv --bar-minutes 5 \
        --features f_a,f_b,... [--model ridge|logistic] [--out DIR] \
        [--n-shuffle 50] [--skip-lgbm] \
        [--write-artifact PATH --tf 5min --confirm-tfs 15min,1h]

Analysis tooling only — never imported by the live code.
"""
from __future__ import annotations

import argparse
import hashlib
import json
import math
import subprocess
import sys
import warnings
from pathlib import Path
from typing import Dict, List, Optional, Sequence, Tuple

import numpy as np
import pandas as pd

PROJECT_ROOT = Path(__file__).resolve().parents[2]
if str(PROJECT_ROOT) not in sys.path:
    sys.path.insert(0, str(PROJECT_ROOT))

from position_sizing import FEATURE_ORDER, FEATURE_VERSION  # noqa: E402  (read-only import)

# ── Pre-registered constants (flag any change loudly) ──────────────────
WINSOR_LO, WINSOR_HI = -1.5, 4.0
LAMBDA_GRID: Tuple[float, ...] = (1.0, 10.0, 100.0, 1000.0)
CV_FOLDS = 5
MIN_TRAIN = 150
MAP_Q_LO, MAP_Q_HI = 0.10, 0.90
MAP_OUT_LO, MAP_OUT_HI = 0.5, 1.5
MAP_CLIP_LO, MAP_CLIP_HI = 0.5, 1.5
MAP_ROUND_TO = 0.05
MEAN_MULT_TOL = 0.05
DD_TOL = 1.10
PF_FLOOR = 1.3
NEWTON_MAX_ITER = 50
NEWTON_TOL = 1e-8

#: The exhaustion-block ablation's fixed feature list.
EXHAUSTION_FEATURES: Tuple[str, ...] = (
    "f_swing_dist_atr",
    "f_ema21_dist_atr",
    "f_ema200_dist_atr",
    "f_exhaustion_count",
    "f_penalty_total",
    "f_tag_bb_overextended",
    "f_tag_rsi_extended",
    "f_tag_no_pullback",
    "f_tag_exhaustion",
)
ISBUY_FEATURES: Tuple[str, ...] = ("f_is_buy",)

#: lightgbm never ships; running the full shuffle null through 200-tree fits
#: costs minutes for a number nobody acts on, so its null is capped here.
#: DEVIATION FROM THE PROTOCOL CONSTANT — printed in the summary when used.
LGBM_MAX_SHUFFLE = 10

LGBM_PARAMS = dict(
    n_estimators=200,
    learning_rate=0.03,
    num_leaves=7,
    min_child_samples=30,
    subsample=0.8,
    colsample_bytree=0.8,
    reg_lambda=10,
)

_LINES: List[str] = []


def say(msg: str = "") -> None:
    """Print and tee into the summary buffer."""
    print(msg)
    _LINES.append(msg)


# ══════════════════════════════════════════════════════════════════════
# Book loading
# ══════════════════════════════════════════════════════════════════════
def load_book(path: Path) -> pd.DataFrame:
    """Read a feature book, parse ``entry_time`` UTC, drop NaN-R rows."""
    df = pd.read_csv(path)
    if "entry_time" not in df.columns or "r_multiple" not in df.columns:
        raise SystemExit(f"{path}: not a trade book (no entry_time/r_multiple)")
    df["entry_time"] = pd.to_datetime(df["entry_time"], utc=True, format="mixed")
    n0 = len(df)
    bad = ~np.isfinite(pd.to_numeric(df["r_multiple"], errors="coerce").to_numpy(float))
    if bad.any():
        say(f"WARNING: dropping {int(bad.sum())}/{n0} rows with non-finite r_multiple")
        df = df.loc[~bad].copy()
    df["r_multiple"] = df["r_multiple"].astype(float)
    df = df.sort_values("entry_time", kind="mergesort").reset_index(drop=True)
    return df


def resolve_features(df: pd.DataFrame, names: Sequence[str], label: str) -> List[str]:
    """Validate a requested feature list against the book and FEATURE_ORDER."""
    missing = [n for n in names if n not in df.columns]
    if missing:
        raise SystemExit(f"{label}: columns absent from book: {', '.join(missing)}")
    unknown = [n for n in names if n.startswith("f_") and n not in FEATURE_ORDER]
    if unknown:
        say(f"NOTE: {label}: {', '.join(unknown)} not in position_sizing.FEATURE_ORDER")
    return list(names)


def sha256_file(path: Path) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        for chunk in iter(lambda: fh.read(1 << 20), b""):
            h.update(chunk)
    return h.hexdigest()


def git_commit() -> str:
    try:
        return subprocess.check_output(
            ["git", "rev-parse", "HEAD"], cwd=PROJECT_ROOT, text=True).strip()
    except Exception:  # pragma: no cover - research tool
        return "unknown"


# ══════════════════════════════════════════════════════════════════════
# Walk-forward geometry
# ══════════════════════════════════════════════════════════════════════
def quarter_starts(first: pd.Timestamp, last: pd.Timestamp) -> List[pd.Timestamp]:
    """Every Jan/Apr/Jul/Oct 1 00:00 UTC in ``[first, last]`` inclusive."""
    out: List[pd.Timestamp] = []
    year = first.year
    while True:
        for month in (1, 4, 7, 10):
            ts = pd.Timestamp(year=year, month=month, day=1, tz="UTC")
            if ts < first:
                continue
            if ts > last:
                return out
            out.append(ts)
        year += 1


def quarter_label(ts: pd.Timestamp) -> str:
    return f"{ts.year}Q{(ts.month - 1) // 3 + 1}"


def build_folds(times: np.ndarray, n_folds: int) -> List[np.ndarray]:
    """Contiguous-by-time index blocks (rows are already time-sorted)."""
    n = len(times)
    edges = np.linspace(0, n, n_folds + 1).astype(int)
    return [np.arange(edges[i], edges[i + 1]) for i in range(n_folds)
            if edges[i + 1] > edges[i]]


# ══════════════════════════════════════════════════════════════════════
# Standardisation
# ══════════════════════════════════════════════════════════════════════
def fit_standardiser(X: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return (mean, scale, keep) from a training matrix.

    NaNs are ignored when computing the mean (they are later filled with it).
    ``scale`` is the population std with zeros replaced by 1.0; ``keep`` is
    False for zero-variance-in-train columns, which are dropped from the fit
    and carry a coefficient of 0.
    """
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        mean = np.nanmean(X, axis=0)
        std = np.nanstd(X, axis=0)
    mean = np.where(np.isfinite(mean), mean, 0.0)
    std = np.where(np.isfinite(std), std, 0.0)
    keep = std > 0
    scale = np.where(keep, std, 1.0)
    return mean, scale, keep


def apply_standardiser(X: np.ndarray, mean: np.ndarray, scale: np.ndarray) -> np.ndarray:
    """Fill NaN with the TRAIN mean, then standardise."""
    Z = np.array(X, dtype=float, copy=True)
    bad = ~np.isfinite(Z)
    if bad.any():
        Z[bad] = np.take(mean, np.where(bad)[1])
    return (Z - mean) / scale


# ══════════════════════════════════════════════════════════════════════
# Models (numpy only)
# ══════════════════════════════════════════════════════════════════════
def _solve(A: np.ndarray, b: np.ndarray) -> np.ndarray:
    try:
        return np.linalg.solve(A, b)
    except np.linalg.LinAlgError:  # pragma: no cover - degenerate design
        return np.linalg.lstsq(A, b, rcond=None)[0]


def fit_ridge(X: np.ndarray, y: np.ndarray, lam: float) -> Tuple[np.ndarray, float]:
    """Closed-form ridge with an UNPENALISED intercept.

    Centring X and y is exactly equivalent to appending a 1-column carrying
    zero penalty. Returns ``(coef, intercept)`` with ``pred = intercept + X @ coef``.
    """
    if X.shape[1] == 0:
        return np.zeros(0), float(y.mean()) if len(y) else 0.0
    xm = X.mean(axis=0)
    ym = float(y.mean())
    Xc = X - xm
    w = _solve(Xc.T @ Xc + lam * np.eye(X.shape[1]), Xc.T @ (y - ym))
    return w, ym - float(xm @ w)


def fit_logistic(X: np.ndarray, y: np.ndarray, lam: float,
                 max_iter: int = NEWTON_MAX_ITER,
                 tol: float = NEWTON_TOL) -> Tuple[np.ndarray, float, float]:
    """L2 logistic via Newton-Raphson, intercept unpenalised.

    Returns ``(coef, intercept, grad_norm)`` where ``grad_norm`` is the
    infinity-norm of the penalised gradient at the returned point.
    """
    n, p = X.shape
    A = np.hstack([np.ones((n, 1)), X])
    pen = np.zeros(p + 1)
    pen[1:] = lam
    beta = np.zeros(p + 1)
    # Warm-start the intercept at the empirical log-odds.
    frac = float(np.clip(y.mean(), 1e-6, 1 - 1e-6))
    beta[0] = math.log(frac / (1 - frac))
    grad = np.zeros(p + 1)
    for _ in range(max_iter):
        eta = A @ beta
        prob = 1.0 / (1.0 + np.exp(-np.clip(eta, -35, 35)))
        grad = A.T @ (y - prob) - pen * beta
        W = np.clip(prob * (1 - prob), 1e-10, None)
        H = A.T @ (A * W[:, None]) + np.diag(pen)
        step = _solve(H, grad)
        beta = beta + step
        if np.max(np.abs(step)) < tol:
            break
    eta = A @ beta
    prob = 1.0 / (1.0 + np.exp(-np.clip(eta, -35, 35)))
    grad = A.T @ (y - prob) - pen * beta
    return beta[1:], float(beta[0]), float(np.max(np.abs(grad)))


def predict_linear(X: np.ndarray, coef: np.ndarray, intercept: float) -> np.ndarray:
    if X.shape[1] == 0:
        return np.full(len(X), intercept, dtype=float)
    return intercept + X @ coef


def sigmoid(z: np.ndarray) -> np.ndarray:
    return 1.0 / (1.0 + np.exp(-np.clip(z, -35, 35)))


def log_loss(y: np.ndarray, p: np.ndarray) -> float:
    p = np.clip(p, 1e-12, 1 - 1e-12)
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))


def choose_lambda(X: np.ndarray, y: np.ndarray, kind: str,
                  grid: Sequence[float] = LAMBDA_GRID,
                  n_folds: int = CV_FOLDS) -> Tuple[float, Dict[float, float]]:
    """5-fold time-contiguous CV inside the training window."""
    folds = build_folds(np.arange(len(y)), n_folds)
    if X.shape[1] == 0 or len(folds) < 2:
        return float(grid[0]), {}
    scores: Dict[float, float] = {}
    for lam in grid:
        vals: List[float] = []
        for hold in folds:
            mask = np.ones(len(y), dtype=bool)
            mask[hold] = False
            if mask.sum() < 2:
                continue
            if kind == "logistic":
                ytr = y[mask]
                if ytr.min() == ytr.max():          # single-class fold
                    p = np.full(len(hold), float(np.clip(ytr.mean(), 1e-6, 1 - 1e-6)))
                else:
                    c, b, _ = fit_logistic(X[mask], ytr, lam)
                    p = sigmoid(predict_linear(X[hold], c, b))
                vals.append(log_loss(y[hold], p))
            else:
                c, b = fit_ridge(X[mask], y[mask], lam)
                resid = y[hold] - predict_linear(X[hold], c, b)
                vals.append(float(np.mean(resid ** 2)))
        scores[float(lam)] = float(np.mean(vals)) if vals else float("inf")
    best = min(scores, key=lambda k: scores[k])
    return best, scores


def lgbm_available() -> bool:
    try:
        import lightgbm  # noqa: F401
        return True
    except Exception:
        return False


def fit_lgbm(X: np.ndarray, y: np.ndarray, seed: int):
    import lightgbm as lgb
    model = lgb.LGBMRegressor(random_state=seed, verbose=-1, **LGBM_PARAMS)
    model.fit(X, y)
    return model


# ══════════════════════════════════════════════════════════════════════
# Size map
# ══════════════════════════════════════════════════════════════════════
def round_to(x, step: float = MAP_ROUND_TO) -> np.ndarray:
    """Round to the nearest ``step`` and kill float dust (0.05 → 2 dp)."""
    decimals = max(0, int(round(-math.log10(step)))) + 1
    return np.round(np.round(np.asarray(x, dtype=float) / step) * step, decimals)


def calibrate_map(z_train: np.ndarray) -> Dict[str, float]:
    """lo/hi/normaliser from TRAIN predictions only."""
    lo = float(np.quantile(z_train, MAP_Q_LO))
    hi = float(np.quantile(z_train, MAP_Q_HI))
    if not np.isfinite(lo) or not np.isfinite(hi) or hi <= lo:
        return {"lo": lo, "hi": hi, "normalizer": 1.0, "degenerate": True}
    raw = MAP_OUT_LO + (MAP_OUT_HI - MAP_OUT_LO) * np.clip((z_train - lo) / (hi - lo), 0, 1)
    mean_raw = float(np.mean(raw))
    norm = 1.0 / mean_raw if mean_raw > 0 else 1.0
    return {"lo": lo, "hi": hi, "normalizer": norm, "degenerate": False}


def apply_map(z: np.ndarray, cal: Dict[str, float]) -> np.ndarray:
    z = np.asarray(z, dtype=float)
    if cal.get("degenerate"):
        return np.ones(len(z))
    lo, hi, norm = cal["lo"], cal["hi"], cal["normalizer"]
    raw = MAP_OUT_LO + (MAP_OUT_HI - MAP_OUT_LO) * np.clip((z - lo) / (hi - lo), 0, 1)
    return round_to(np.clip(raw * norm, MAP_CLIP_LO, MAP_CLIP_HI))


# ══════════════════════════════════════════════════════════════════════
# Metrics
# ══════════════════════════════════════════════════════════════════════
def max_drawdown(r: np.ndarray) -> float:
    """0-prepended cumulative-R drawdown (gpu_backtest.compute_metrics)."""
    r = np.asarray(r, dtype=float)
    if len(r) == 0:
        return 0.0
    cum = np.concatenate(([0.0], np.cumsum(r)))
    return float((np.maximum.accumulate(cum) - cum).max())


def profit_factor(r: np.ndarray) -> float:
    r = np.asarray(r, dtype=float)
    gp = float(r[r > 0].sum()) if (r > 0).any() else 0.0
    gl = float(-r[r < 0].sum()) if (r < 0).any() else 0.0
    return gp / gl if gl > 0 else float("inf")


def multiplier_hist(m: np.ndarray) -> pd.DataFrame:
    """Counts in 0.1-wide bins from 0.5 to 1.5.

    ``np.histogram`` already closes the final bin on the right, so an
    exactly-1.5 multiplier lands in [1.4, 1.5] once and the counts sum to
    ``len(m)``.
    """
    edges = np.round(np.arange(MAP_OUT_LO, MAP_OUT_HI + 1e-9, 0.1), 2)
    counts, _ = np.histogram(np.asarray(m, dtype=float), bins=edges)
    counts = counts.astype(int)
    return pd.DataFrame({
        "bin_lo": edges[:-1], "bin_hi": edges[1:], "count": counts,
        "pct": np.where(len(m), counts / max(len(m), 1) * 100, 0.0),
    })


# ══════════════════════════════════════════════════════════════════════
# Walk-forward driver
# ══════════════════════════════════════════════════════════════════════
def walk_forward(df: pd.DataFrame, features: Sequence[str], model: str,
                 bar_minutes: int, burn_in_months: int, purge_bars: int,
                 seed: int, y_train_override: Optional[np.ndarray] = None,
                 verbose: bool = True) -> Dict:
    """One full expanding walk-forward pass.

    ``y_train_override`` (full-book length, aligned to ``df``) replaces the
    TRAINING labels only — this is how the shuffle null is run. Test labels
    always come from the real ``r_multiple``.
    """
    times = df["entry_time"].to_numpy()
    r_all = df["r_multiple"].to_numpy(float)
    r_train_src = r_all if y_train_override is None else np.asarray(y_train_override, float)
    X_all = df[list(features)].apply(pd.to_numeric, errors="coerce").to_numpy(float) \
        if features else np.zeros((len(df), 0))

    book_start = df["entry_time"].iloc[0]
    book_end = df["entry_time"].iloc[-1]
    burn_end = book_start + pd.DateOffset(months=burn_in_months)
    purge = pd.Timedelta(minutes=purge_bars * bar_minutes)

    qs = [q for q in quarter_starts(book_start, book_end) if q >= burn_end]
    rows: List[Dict] = []
    per_q: List[Dict] = []
    coef_rows: List[Dict] = []
    skipped: List[str] = []
    dropped_seen: set = set()

    for q_start in qs:
        q_end = q_start + pd.DateOffset(months=3)
        train_mask = df["entry_time"] < (q_start - purge)
        test_mask = (df["entry_time"] >= q_start) & (df["entry_time"] < q_end)
        n_tr, n_te = int(train_mask.sum()), int(test_mask.sum())
        if n_tr < MIN_TRAIN or n_te == 0:
            skipped.append(f"{quarter_label(q_start)} (train={n_tr}, test={n_te})")
            continue

        tr = np.flatnonzero(train_mask.to_numpy())
        te = np.flatnonzero(test_mask.to_numpy())
        Xtr_raw, Xte_raw = X_all[tr], X_all[te]
        mean, scale, keep = fit_standardiser(Xtr_raw)
        dropped = [features[i] for i in np.flatnonzero(~keep)]
        if verbose and dropped:
            tag = (quarter_label(q_start), tuple(dropped))
            if tag not in dropped_seen:
                dropped_seen.add(tag)
                say(f"  WARNING {quarter_label(q_start)}: zero variance in train, "
                    f"dropped from fit (coef 0): {', '.join(dropped)}")
        Xtr = apply_standardiser(Xtr_raw, mean, scale)[:, keep]
        Xte = apply_standardiser(Xte_raw, mean, scale)[:, keep]

        r_tr = r_train_src[tr]
        lam: Optional[float] = None
        if model == "logistic":
            y_tr = (r_tr > 0).astype(float)
            lam, _ = choose_lambda(Xtr, y_tr, "logistic")
            coef, b, _ = fit_logistic(Xtr, y_tr, lam)
            z_tr = sigmoid(predict_linear(Xtr, coef, b))
            z_te = sigmoid(predict_linear(Xte, coef, b))
        elif model == "lgbm":
            y_tr = np.clip(r_tr, WINSOR_LO, WINSOR_HI)
            gbm = fit_lgbm(Xtr, y_tr, seed)
            z_tr = np.asarray(gbm.predict(Xtr), dtype=float)
            z_te = np.asarray(gbm.predict(Xte), dtype=float)
            coef, b = np.zeros(int(keep.sum())), 0.0
        else:  # ridge / isbuy / exhaustion
            y_tr = np.clip(r_tr, WINSOR_LO, WINSOR_HI)
            lam, _ = choose_lambda(Xtr, y_tr, "ridge")
            coef, b = fit_ridge(Xtr, y_tr, lam)
            z_tr = predict_linear(Xtr, coef, b)
            z_te = predict_linear(Xte, coef, b)

        cal = calibrate_map(z_tr)
        m_te = apply_map(z_te, cal)

        if model != "lgbm":
            full_coef = np.zeros(len(features))
            full_coef[np.flatnonzero(keep)] = coef
            for name, c in zip(features, full_coef):
                coef_rows.append({"quarter": quarter_label(q_start), "feature": name,
                                  "coef_std": float(c), "lambda": lam,
                                  "dropped": name in dropped})
            coef_rows.append({"quarter": quarter_label(q_start), "feature": "_intercept",
                              "coef_std": float(b), "lambda": lam, "dropped": False})

        r_te = r_all[te]
        for i, idx in enumerate(te):
            rows.append({
                "entry_time": df["entry_time"].iloc[idx],
                "r_multiple": float(r_te[i]),
                "mult": float(m_te[i]),
                "r_sized": float(r_te[i] * m_te[i]),
                "quarter": quarter_label(q_start),
                "z": float(z_te[i]),
            })
        per_q.append({
            "quarter": quarter_label(q_start),
            "q_start": q_start, "n_train": n_tr, "n_test": n_te,
            "flat_r": float(r_te.sum()), "sized_r": float((r_te * m_te).sum()),
            "mean_m": float(m_te.mean()), "lambda": lam,
            "map_lo": cal["lo"], "map_hi": cal["hi"], "normalizer": cal["normalizer"],
        })

    oos = pd.DataFrame(rows)
    return {"oos": oos, "quarters": pd.DataFrame(per_q),
            "coefs": pd.DataFrame(coef_rows), "skipped": skipped,
            "quarter_starts": qs}


def evaluate(oos: pd.DataFrame) -> Dict:
    """Headline OOS numbers for one model's walk-forward output."""
    if oos.empty:
        return {"n": 0, "flat_r": 0.0, "sized_r": 0.0, "lift": 0.0}
    r = oos["r_multiple"].to_numpy(float)
    rs = oos["r_sized"].to_numpy(float)
    m = oos["mult"].to_numpy(float)
    med = oos["entry_time"].median()
    first = (oos["entry_time"] < med).to_numpy()
    second = ~first
    return {
        "n": int(len(oos)),
        "flat_r": float(r.sum()), "sized_r": float(rs.sum()),
        "lift": float(rs.sum() - r.sum()),
        "flat_pf": profit_factor(r), "sized_pf": profit_factor(rs),
        "flat_dd": max_drawdown(r), "sized_dd": max_drawdown(rs),
        "mean_m": float(m.mean()), "min_m": float(m.min()), "max_m": float(m.max()),
        "median_time": med,
        "n_h1": int(first.sum()), "n_h2": int(second.sum()),
        "flat_h1": float(r[first].sum()), "flat_h2": float(r[second].sum()),
        "sized_h1": float(rs[first].sum()), "sized_h2": float(rs[second].sum()),
    }


def acceptance(ev: Dict, real_lift: float, shuffle_p95: float) -> Dict:
    """The six pre-registered clauses."""
    c1 = ev["sized_h1"] > ev["flat_h1"] and ev["sized_h2"] > ev["flat_h2"]
    c2 = ev["sized_r"] > ev["flat_r"]
    c3 = ev["sized_dd"] <= DD_TOL * ev["flat_dd"]
    c4 = (ev["sized_pf"] >= PF_FLOOR) and ev["sized_h1"] > 0 and ev["sized_h2"] > 0
    c5 = real_lift > shuffle_p95 if np.isfinite(shuffle_p95) else False
    c6 = (abs(ev["mean_m"] - 1.0) <= MEAN_MULT_TOL
          and ev["min_m"] >= MAP_CLIP_LO - 1e-9 and ev["max_m"] <= MAP_CLIP_HI + 1e-9)
    clauses = {"1_both_halves": bool(c1), "2_overall": bool(c2), "3_drawdown": bool(c3),
               "4_pf_and_halves": bool(c4), "5_beats_shuffle": bool(c5),
               "6_multiplier_bounds": bool(c6)}
    clauses["overall"] = bool(all(clauses.values()))
    return clauses


def shuffle_null(df: pd.DataFrame, features: Sequence[str], model: str,
                 bar_minutes: int, burn_in_months: int, purge_bars: int,
                 n_shuffle: int, seed: int) -> np.ndarray:
    """Lift distribution under permuted TRAINING labels (test labels real)."""
    rng = np.random.default_rng(seed)
    r_all = df["r_multiple"].to_numpy(float)
    lifts: List[float] = []
    for _ in range(n_shuffle):
        r_perm = r_all[rng.permutation(len(r_all))]
        res = walk_forward(df, features, model, bar_minutes, burn_in_months,
                           purge_bars, seed, y_train_override=r_perm, verbose=False)
        ev = evaluate(res["oos"])
        lifts.append(ev["lift"])
    return np.asarray(lifts, dtype=float)


# ══════════════════════════════════════════════════════════════════════
# Final fit + artifact
# ══════════════════════════════════════════════════════════════════════
def final_fit(df: pd.DataFrame, features: Sequence[str], model: str,
              seed: int) -> Dict:
    """Fit on ALL trades; calibrate the map on in-sample predictions."""
    X_raw = df[list(features)].apply(pd.to_numeric, errors="coerce").to_numpy(float)
    r = df["r_multiple"].to_numpy(float)
    mean, scale, keep = fit_standardiser(X_raw)
    X = apply_standardiser(X_raw, mean, scale)[:, keep]
    if model == "logistic":
        y = (r > 0).astype(float)
        lam, _ = choose_lambda(X, y, "logistic")
        coef, b, gnorm = fit_logistic(X, y, lam)
        z = sigmoid(predict_linear(X, coef, b))
        target = "r_positive"
    else:
        y = np.clip(r, WINSOR_LO, WINSOR_HI)
        lam, _ = choose_lambda(X, y, "ridge")
        coef, b = fit_ridge(X, y, lam)
        gnorm = float("nan")
        z = predict_linear(X, coef, b)
        target = "r_winsorised"
    full_coef = np.zeros(len(features))
    full_coef[np.flatnonzero(keep)] = coef
    cal = calibrate_map(z)
    return {"mean": mean, "scale": scale, "keep": keep, "coef": full_coef,
            "intercept": float(b), "lambda": float(lam), "z": z,
            "mult": apply_map(z, cal), "map": cal, "target": target,
            "grad_norm": gnorm}


def predict_multiplier(feature_dict: Dict[str, float], artifact: Dict) -> float:
    """Runtime formula: standardise → linear score → size map → 0.05 grid."""
    names = artifact["features"]
    mean = artifact["mean"]
    scale = artifact["scale"]
    coef = artifact["coef"]
    z = float(artifact["intercept"])
    for i, name in enumerate(names):
        raw = feature_dict.get(name, mean[i])
        if raw is None or not np.isfinite(float(raw)):
            raw = mean[i]
        z += float(coef[i]) * (float(raw) - float(mean[i])) / float(scale[i])
    if artifact["model"] == "logistic":
        z = 1.0 / (1.0 + math.exp(-max(-35.0, min(35.0, z))))
    mp = artifact["map"]
    if mp["hi"] <= mp["lo"]:
        return 1.0
    frac = min(1.0, max(0.0, (z - mp["lo"]) / (mp["hi"] - mp["lo"])))
    raw_m = mp["out_lo"] + (mp["out_hi"] - mp["out_lo"]) * frac
    val = min(mp["clip_hi"], max(mp["clip_lo"], raw_m * mp["normalizer"]))
    step = mp["round_to"]
    return float(round(round(val / step) * step, 2))


def build_artifact(fit: Dict, features: Sequence[str], model: str, df: pd.DataFrame,
                   ev: Dict, clauses: Dict, real_lift: float, shuffle: np.ndarray,
                   n_shuffle: int, book: Path, tf: str, confirm_tfs: List[str],
                   burn_in_months: int, purge_bars: int, spread: float) -> Dict:
    return {
        "schema_version": 1,
        "feature_version": FEATURE_VERSION,
        "pair": "XAU/USD",
        "primary_tf": tf,
        "confirmation_tfs": confirm_tfs,
        "model": model,
        "target": fit["target"],
        "features": list(features),
        "mean": [float(v) for v in fit["mean"]],
        "scale": [float(v) for v in fit["scale"]],
        "coef": [float(v) for v in fit["coef"]],
        "intercept": float(fit["intercept"]),
        "map": {
            "kind": "clipped_linear_on_quantiles",
            "lo": float(fit["map"]["lo"]), "hi": float(fit["map"]["hi"]),
            "out_lo": MAP_OUT_LO, "out_hi": MAP_OUT_HI,
            "normalizer": float(fit["map"]["normalizer"]),
            "clip_lo": MAP_CLIP_LO, "clip_hi": MAP_CLIP_HI,
            "round_to": MAP_ROUND_TO,
        },
        "training": {
            "window": [df["entry_time"].iloc[0].isoformat(),
                       df["entry_time"].iloc[-1].isoformat()],
            "spread": spread,
            "n_train_trades": int(len(df)),
            "n_oos_trades": int(ev.get("n", 0)),
            "burn_in_months": burn_in_months,
            "purge_bars": purge_bars,
            "refit": "calendar_quarter",
            "final_fit": "full_window",
            "lambda": float(fit["lambda"]),
            "git_commit": git_commit(),
            "book_sha256": sha256_file(book),
        },
        "acceptance": {
            **{k: v for k, v in clauses.items()},
            "flat_r": ev.get("flat_r"), "sized_r": ev.get("sized_r"),
            "sized_pf": ev.get("sized_pf"), "flat_dd": ev.get("flat_dd"),
            "sized_dd": ev.get("sized_dd"),
            "sized_h1": ev.get("sized_h1"), "sized_h2": ev.get("sized_h2"),
            "flat_h1": ev.get("flat_h1"), "flat_h2": ev.get("flat_h2"),
            "real_lift": float(real_lift),
            "shuffle_p95_lift": float(np.percentile(shuffle, 95)) if len(shuffle) else None,
            "n_shuffle": int(n_shuffle),
            "mean_multiplier": ev.get("mean_m"),
        },
    }


# ══════════════════════════════════════════════════════════════════════
# Reporting
# ══════════════════════════════════════════════════════════════════════
def fmt(x, nd=2, width=0) -> str:
    if x is None:
        return "n/a".rjust(width)
    if isinstance(x, float) and not np.isfinite(x):
        return ("inf" if x > 0 else "-inf").rjust(width)
    s = f"{x:.{nd}f}" if isinstance(x, (int, float, np.floating)) else str(x)
    return s.rjust(width) if width else s


def print_quarters(q: pd.DataFrame) -> None:
    say(f"{'quarter':<9}{'n_tr':>6}{'n_te':>6}{'flat R':>10}{'sized R':>10}"
        f"{'mean m':>9}{'lambda':>9}")
    for _, row in q.iterrows():
        say(f"{row['quarter']:<9}{int(row['n_train']):>6}{int(row['n_test']):>6}"
            f"{fmt(row['flat_r'], 2, 10)}{fmt(row['sized_r'], 2, 10)}"
            f"{fmt(row['mean_m'], 3, 9)}"
            f"{('n/a' if row['lambda'] is None else fmt(row['lambda'], 0)):>9}")


def print_hist(m: np.ndarray) -> None:
    h = multiplier_hist(m)
    say(f"{'bin':<12}{'count':>8}{'pct':>8}")
    for _, row in h.iterrows():
        bar = "#" * int(round(row["pct"] / 2))
        say(f"{row['bin_lo']:.1f}–{row['bin_hi']:.1f}{'':<4}{int(row['count']):>8}"
            f"{row['pct']:>7.1f}%  {bar}")


def print_coef_signs(coefs: pd.DataFrame, label: str) -> None:
    """Mean standardised coefficient + sign stability across quarters."""
    if coefs.empty:
        return
    d = coefs[coefs["feature"] != "_intercept"]
    say(f"  standardised coefficients across quarters ({label}):")
    say(f"    {'feature':<24}{'mean':>10}{'min':>10}{'max':>10}{'+q':>5}{'-q':>5}  sign-stable")
    for name, grp in d.groupby("feature", sort=False):
        c = grp["coef_std"].to_numpy(float)
        pos, neg = int((c > 0).sum()), int((c < 0).sum())
        stable = "yes" if (pos == 0 or neg == 0) and (pos + neg) > 0 else "no"
        say(f"    {name:<24}{c.mean():>10.4f}{c.min():>10.4f}{c.max():>10.4f}"
            f"{pos:>5}{neg:>5}  {stable}")


def print_comparison(table: List[Dict]) -> None:
    say(f"{'model':<12}{'n':>6}{'flat R':>10}{'sized R':>10}{'lift':>9}"
        f"{'PF':>7}{'DD':>8}{'mean m':>8}{'p95 null':>10}{'rank%':>7}  pass")
    for row in table:
        say(f"{row['model']:<12}{row['n']:>6}{fmt(row['flat_r'], 2, 10)}"
            f"{fmt(row['sized_r'], 2, 10)}{fmt(row['lift'], 2, 9)}"
            f"{fmt(row['sized_pf'], 2, 7)}{fmt(row['sized_dd'], 2, 8)}"
            f"{fmt(row['mean_m'], 3, 8)}{fmt(row['shuffle_p95'], 2, 10)}"
            f"{fmt(row['rank_pct'], 0, 7)}  {'PASS' if row['pass'] else 'FAIL'}")


# ══════════════════════════════════════════════════════════════════════
# main
# ══════════════════════════════════════════════════════════════════════
def main() -> None:
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--book", required=True, type=Path)
    ap.add_argument("--bar-minutes", required=True, type=int)
    ap.add_argument("--features", required=True,
                    help="comma-separated f_* column names for the main model")
    ap.add_argument("--model", default="ridge", choices=("ridge", "logistic"))
    ap.add_argument("--burn-in-months", type=int, default=12)
    ap.add_argument("--purge-bars", type=int, default=300)
    ap.add_argument("--n-shuffle", type=int, default=50)
    ap.add_argument("--seed", type=int, default=20260905)
    ap.add_argument("--skip-lgbm", action="store_true")
    ap.add_argument("--out", type=Path, default=None)
    ap.add_argument("--write-artifact", type=Path, default=None)
    ap.add_argument("--tf", default=None, help="primary TF label for the artifact")
    ap.add_argument("--confirm-tfs", default="15min,1h")
    ap.add_argument("--spread", type=float, default=0.65,
                    help="cost the book was produced at (recorded in the artifact)")
    args = ap.parse_args()

    if args.write_artifact and not args.tf:
        raise SystemExit("--write-artifact requires --tf")

    say("=" * 96)
    say("Stage-D stake-sizing walk-forward — RESEARCH ONLY")
    say("=" * 96)
    say(f"book            : {args.book}")
    df = load_book(args.book)
    say(f"trades          : {len(df)}  "
        f"({df['entry_time'].iloc[0]} → {df['entry_time'].iloc[-1]})")
    say(f"bar minutes     : {args.bar_minutes}   purge = {args.purge_bars} bars = "
        f"{args.purge_bars * args.bar_minutes} min "
        f"({args.purge_bars * args.bar_minutes / 1440:.2f} d)")
    say(f"burn-in         : {args.burn_in_months} months  "
        f"(train-only until {df['entry_time'].iloc[0] + pd.DateOffset(months=args.burn_in_months)})")
    say(f"model / seed    : {args.model} / {args.seed}   shuffles: {args.n_shuffle}")
    say(f"winsor / lambda : [{WINSOR_LO}, {WINSOR_HI}] / {list(LAMBDA_GRID)}  "
        f"CV folds {CV_FOLDS} (time-contiguous)")
    say(f"map             : m = round_{MAP_ROUND_TO}(clip(norm × (0.5 + clip((z−lo)/(hi−lo),0,1)"
        f"), {MAP_CLIP_LO}, {MAP_CLIP_HI}))  lo/hi = train q{MAP_Q_LO}/q{MAP_Q_HI}")
    say("")

    main_feats = resolve_features(df, [f.strip() for f in args.features.split(",") if f.strip()],
                                  "--features")
    say(f"main features ({len(main_feats)}): {', '.join(main_feats)}")

    specs: List[Tuple[str, List[str], str]] = [(args.model, main_feats, args.model)]
    for label, feats in (("isbuy", ISBUY_FEATURES), ("exhaustion", EXHAUSTION_FEATURES)):
        have = [f for f in feats if f in df.columns]
        if len(have) != len(feats):
            say(f"NOTE: {label} baseline missing "
                f"{', '.join(f for f in feats if f not in df.columns)} — using {len(have)}")
        if have:
            specs.append((label, have, "ridge"))
    use_lgbm = (not args.skip_lgbm) and lgbm_available()
    if not args.skip_lgbm and not use_lgbm:
        say("lightgbm not installed — upper bound skipped")
    if use_lgbm:
        all_f = [c for c in df.columns if c.startswith("f_")]
        specs.append(("lgbm", all_f, "lgbm"))
        say(f"lightgbm upper bound: {len(all_f)} f_* columns")
    say("")

    results: Dict[str, Dict] = {}
    table: List[Dict] = []
    for label, feats, kind in specs:
        say("-" * 96)
        say(f"MODEL: {label}   ({kind}, {len(feats)} feature(s))")
        say("-" * 96)
        res = walk_forward(df, feats, kind, args.bar_minutes, args.burn_in_months,
                           args.purge_bars, args.seed, verbose=True)
        if res["skipped"]:
            say(f"  skipped quarters: {'; '.join(res['skipped'])}")
        if res["oos"].empty:
            say("  no OOS quarters — nothing to evaluate")
            continue
        ev = evaluate(res["oos"])
        n_sh = args.n_shuffle if kind != "lgbm" else min(args.n_shuffle, LGBM_MAX_SHUFFLE)
        if kind == "lgbm" and n_sh != args.n_shuffle:
            say(f"  NOTE (protocol deviation, flagged): lgbm shuffle null capped at "
                f"{n_sh} repeats (never ships)")
        sh = shuffle_null(df, feats, kind, args.bar_minutes, args.burn_in_months,
                          args.purge_bars, n_sh, args.seed) if n_sh > 0 else np.array([])
        p95 = float(np.percentile(sh, 95)) if len(sh) else float("nan")
        rank = float((sh < ev["lift"]).mean() * 100) if len(sh) else float("nan")
        cl = acceptance(ev, ev["lift"], p95)

        say(f"  OOS quarters    : {len(res['quarters'])} "
            f"({res['quarters']['quarter'].iloc[0]} → {res['quarters']['quarter'].iloc[-1]}), "
            f"{ev['n']} trades")
        say(f"  flat  R         : {ev['flat_r']:+.2f}   PF {fmt(ev['flat_pf'])}   "
            f"maxDD {ev['flat_dd']:.2f}")
        say(f"  sized R         : {ev['sized_r']:+.2f}   PF {fmt(ev['sized_pf'])}   "
            f"maxDD {ev['sized_dd']:.2f}")
        say(f"  lift            : {ev['lift']:+.2f} R")
        say(f"  halves (split {ev['median_time']}): "
            f"h1 flat {ev['flat_h1']:+.2f} / sized {ev['sized_h1']:+.2f} (n={ev['n_h1']}) | "
            f"h2 flat {ev['flat_h2']:+.2f} / sized {ev['sized_h2']:+.2f} (n={ev['n_h2']})")
        say(f"  multiplier      : mean {ev['mean_m']:.3f}  min {ev['min_m']:.2f}  "
            f"max {ev['max_m']:.2f}")
        say(f"  shuffle null    : n={len(sh)}  mean {fmt(float(sh.mean()) if len(sh) else None)}  "
            f"p95 {fmt(p95)}  real-lift percentile rank {fmt(rank, 0)}%")
        say("  acceptance:")
        for k, v in cl.items():
            if k == "overall":
                continue
            say(f"    {k:<22} {'PASS' if v else 'FAIL'}")
        say(f"    {'OVERALL':<22} {'PASS' if cl['overall'] else 'FAIL'}")
        say("  per-quarter:")
        print_quarters(res["quarters"])
        print_coef_signs(res["coefs"], label)
        say("")

        results[label] = {"res": res, "ev": ev, "shuffle": sh, "p95": p95,
                          "rank": rank, "clauses": cl, "features": feats, "kind": kind}
        table.append({"model": label, "n": ev["n"], "flat_r": ev["flat_r"],
                      "sized_r": ev["sized_r"], "lift": ev["lift"],
                      "sized_pf": ev["sized_pf"], "sized_dd": ev["sized_dd"],
                      "mean_m": ev["mean_m"], "shuffle_p95": p95, "rank_pct": rank,
                      "pass": cl["overall"]})

    if args.model not in results:
        raise SystemExit("the main model produced no OOS quarters — check burn-in / "
                         "book window / MIN_TRAIN")

    say("=" * 96)
    say("COMPARISON")
    say("=" * 96)
    print_comparison(table)
    say("")
    say(f"OOS multiplier histogram ({args.model}):")
    print_hist(results[args.model]["res"]["oos"]["mult"].to_numpy())
    say("")

    main_label = args.model
    beats = None
    if main_label in results and "isbuy" in results:
        beats = results[main_label]["ev"]["lift"] > results["isbuy"]["ev"]["lift"]
        say(f"beats isbuy baseline: {'yes' if beats else 'no'}  "
            f"({results[main_label]['ev']['lift']:+.2f} R vs "
            f"{results['isbuy']['ev']['lift']:+.2f} R)")
    say(f"OVERALL ACCEPTANCE ({main_label}): "
        f"{'PASS' if results[main_label]['clauses']['overall'] else 'FAIL'}")
    say("")

    # ── artifact ──────────────────────────────────────────────────────
    artifact = None
    if args.write_artifact:
        say("-" * 96)
        say("FINAL FIT (all trades) + ARTIFACT")
        say("-" * 96)
        fit = final_fit(df, main_feats, args.model, args.seed)
        say(f"  lambda (CV on full window): {fit['lambda']:g}")
        if np.isfinite(fit["grad_norm"]):
            say(f"  logistic |grad|_inf at convergence: {fit['grad_norm']:.3e}")
        say(f"  map lo/hi/normaliser: {fit['map']['lo']:.5f} / {fit['map']['hi']:.5f} / "
            f"{fit['map']['normalizer']:.5f}")
        artifact = build_artifact(
            fit, main_feats, args.model, df, results[main_label]["ev"],
            results[main_label]["clauses"], results[main_label]["ev"]["lift"],
            results[main_label]["shuffle"], args.n_shuffle, args.book, args.tf,
            [t.strip() for t in args.confirm_tfs.split(",") if t.strip()],
            args.burn_in_months, args.purge_bars, args.spread)
        # Self-check: the runtime formula must reproduce the in-sample map exactly.
        rows = df[main_feats].to_dict("records")
        repro = np.array([predict_multiplier(r, artifact) for r in rows])
        exact = bool(np.array_equal(repro, fit["mult"]))
        say(f"  predict_multiplier reproduces in-sample multipliers exactly: "
            f"{'YES' if exact else 'NO'}  (max |diff| "
            f"{float(np.max(np.abs(repro - fit['mult']))):.3e})")
        if not exact:
            raise SystemExit("predict_multiplier does not reproduce the final fit — abort")
        args.write_artifact.parent.mkdir(parents=True, exist_ok=True)
        args.write_artifact.write_text(json.dumps(artifact, indent=2, default=str))
        say(f"  wrote {args.write_artifact}")
        say("")

    # ── outputs ───────────────────────────────────────────────────────
    if args.out:
        out = args.out
        out.mkdir(parents=True, exist_ok=True)
        main_res = results[main_label]
        main_res["res"]["oos"].to_csv(out / "oos_trades.csv", index=False)
        main_res["res"]["quarters"].to_csv(out / "quarters.csv", index=False)
        coefs = pd.concat(
            [r["res"]["coefs"].assign(model=lbl) for lbl, r in results.items()
             if not r["res"]["coefs"].empty],
            ignore_index=True) if results else pd.DataFrame()
        coefs.to_csv(out / "coefs_by_quarter.csv", index=False)
        pd.DataFrame({
            "model": np.repeat([lbl for lbl in results], [len(r["shuffle"]) for r in results.values()]),
            "lift": np.concatenate([r["shuffle"] for r in results.values()])
            if results else np.array([]),
        }).to_csv(out / "shuffle_lifts.csv", index=False)
        multiplier_hist(main_res["res"]["oos"]["mult"].to_numpy()).to_csv(
            out / "multiplier_hist.csv", index=False)
        payload = {
            "book": str(args.book), "book_sha256": sha256_file(args.book),
            "git_commit": git_commit(), "bar_minutes": args.bar_minutes,
            "model": args.model, "features": main_feats, "seed": args.seed,
            "n_shuffle": args.n_shuffle, "burn_in_months": args.burn_in_months,
            "purge_bars": args.purge_bars,
            "beats_isbuy": beats,
            "models": {
                lbl: {
                    "kind": r["kind"], "features": r["features"],
                    "evaluation": {k: (None if isinstance(v, float) and not np.isfinite(v)
                                       else (str(v) if isinstance(v, pd.Timestamp) else v))
                                   for k, v in r["ev"].items()},
                    "shuffle": {"n": int(len(r["shuffle"])),
                                "mean": float(r["shuffle"].mean()) if len(r["shuffle"]) else None,
                                "p95": None if not np.isfinite(r["p95"]) else r["p95"],
                                "real_lift_percentile": None if not np.isfinite(r["rank"])
                                else r["rank"]},
                    "acceptance": r["clauses"],
                } for lbl, r in results.items()
            },
            "artifact": artifact,
        }
        (out / "results.json").write_text(json.dumps(payload, indent=2, default=str))
        (out / "summary.txt").write_text("\n".join(_LINES) + "\n")
        print(f"\nwrote outputs to {out}")


if __name__ == "__main__":
    main()