Wire NNUE eval into search, build the training pipeline, and add a fast

move/eval-agreement benchmark against Stockfish (tests/lop)

Implements the architecture decisions in src/docs/NNUE.md: HalfKAv2_hm
feature set, int16/int8 quantization, per-ply accumulator backup for
king moves, and a lazy-eval-ready integration into search.c/eval.c.
Adds the full annotate_pgns/ training pipeline (extraction, K-fit/
labeling, quiet filtering, tensor building, training, export) plus an
end-to-end orchestrator (run_nnue_pipeline.py) and validation tooling
(sanity_check_net.py, compare_vs_stockfish_net.py). Adds tests/lop:
a 10k-position list-of-positions benchmark comparing engine eval rank
and move choice against pre-computed Stockfish/lc0 references, much
cheaper than a full self-play match for day-to-day progress checks.

First sd10-matched self-play result against classic HCE: LLR +2.28
(H1=+2.94), -72=119+169 over 360 games, confirming the net is
"smarter per node" at matched depth. A real-time (--st 4) match is
currently weaker (candidate loses ~3 plies of depth to a ~3.3-3.7x
nps gap against an unrecalibrated search), consistent with
NNUE.md's adoption-criterion section -- pruning/reduction constants
still need retuning for the new eval's score scale before a fair
real-time comparison is possible.

Co-Authored-By: Claude Sonnet 5 <[email protected]>

18 files changed, 2803 insertions, 144 deletions

diff --git a/annotate_pgns/GNUmakefile b/annotate_pgns/GNUmakefile
new file mode 100644
index 0000000..29fd78d
--- /dev/null
+++ b/annotate_pgns/GNUmakefile
@@ -0,0 +1,54 @@
+# Packages just the NNUE training pipeline's Python driver + its
+# dependencies (NOT the input PGN -- that's passed to
+# run_nnue_pipeline.py's --pgn-glob at runtime, shipped separately) into
+# a single tar.bz2 for copying to a rented GPU box. See docs/NNUE.md.
+#
+# Usage:
+#   make -C annotate_pgns bundle
+#   make -C annotate_pgns clean-bundle
+
+BUNDLE := nnue_pipeline_bundle.tar.bz2
+
+# This directory's own driver + everything it imports.
+PIPELINE_FILES := \
+	run_nnue_pipeline.py \
+	extract_selfplay_positions.py \
+	split_holdout.py \
+	fit_k_and_label.py \
+	quiet_filter_games.py \
+	quiet_filter.py \
+	quiet_common.py \
+	build_training_tensors.py \
+	nnue_feature_spec.py \
+	nnue_model.py \
+	nnue_train.py \
+	export_nnue_weights.py \
+	sample_egtb_positions.py
+
+# quiet_common.py reaches into ../src/eval_tune for its SEE logic
+# (tune_eval_dna.py's static_exchange_eval/has_hanging_material), and
+# dna_diff.py parses ../src/eval.c at import time for the DNA term
+# list -- all three are genuine runtime dependencies, not optional
+# extras, even though they live outside this directory. See
+# quiet_common.py's own docstring.
+EXTERNAL_FILES := \
+	../src/eval_tune/tune_eval_dna.py \
+	../src/eval_tune/dna_diff.py \
+	../src/eval.c
+
+ALL_FILES := $(PIPELINE_FILES) $(EXTERNAL_FILES)
+
+.PHONY: bundle clean-bundle
+
+bundle: $(BUNDLE)
+
+$(BUNDLE): $(ALL_FILES)
+	tar -cjf $@ -C .. \
+		$(addprefix annotate_pgns/,$(PIPELINE_FILES)) \
+		$(patsubst ../%,%,$(EXTERNAL_FILES))
+	@echo "Wrote $@ ($$(du -h $@ | cut -f1)) -- the input PGN is NOT" \
+	      "included; pass it separately via run_nnue_pipeline.py's" \
+	      "--pgn-glob on the rented box."
+
+clean-bundle:
+	rm -f $(BUNDLE)
diff --git a/annotate_pgns/annotate.py b/annotate_pgns/annotate.py
index 240b55f..1f8d02d 100644
--- a/annotate_pgns/annotate.py
+++ b/annotate_pgns/annotate.py
@@ -29,8 +29,13 @@ import chess.pgn
 from quiet_common import position_is_quiet
 
 TELLOTHERS_RE = re.compile(
-    r"^tellothers d(\d+), ([+-]?[\d.]+), [\d.]+s,\s*\d+knps, PV=(.*)$"
+    r"^tellothers d(\d+), ([+-]?[\w.]+), [\d.]+s,\s*\d+knps, PV=(.*)$"
 )
+FLOAT_RE = re.compile(r"^[+-]?\d+\.\d+$")
+MATE_RE = re.compile(r"^([+-])MATE(\d*)$")  # digit(s) optional -- the
+# "FORCED MATE --> stop searching now" immediate-short-circuit path
+# (search.c) prints a bare "+MATE"/"-MATE" with no distance number at
+# all, not just the ordinary "+MATE<N>" case.
 
 
 def load_games(input_paths, max_games=None):
@@ -122,17 +127,52 @@ def run_shard(engine_path, script_path, out_path, hash_size, cpus):
 
 
 def parse_shard_output(out_path):
-    """Return a list of (depth, score, pv_string) in script order."""
+    """Return a list of (depth, score_token, pv_string) in script
+    order. score_token is kept as the engine's raw tellothers token
+    (an ordinary "+1.23"-style float string, or "+MATE<N>"/"-MATE<N>")
+    rather than eagerly parsed to float here -- a MATE token has no
+    pawns magnitude to parse in the first place, and silently dropping
+    any line that doesn't match a plain float is exactly the bug that
+    used to desync this function's output from `index` (every position
+    after a dropped mate/EGTB line in a shard got the wrong score) --
+    see docs/NNUE.md's "Bug found and fixed in annotate.py's output
+    parser" status entry, previously patched only in the separate
+    extract_positions.py reader, not here."""
     results = []
     with open(out_path, encoding="utf-8", errors="replace") as f:
         for line in f:
             m = TELLOTHERS_RE.match(line.rstrip("\n"))
             if m:
-                depth, score, pv = m.groups()
-                results.append((int(depth), float(score), pv.strip()))
+                depth, score_tok, pv = m.groups()
+                results.append((int(depth), score_tok, pv.strip()))
     return results
 
 
+def _format_comment(depth, score_tok, pv):
+    """Render one (depth, score_token, pv) result as a PGN comment in
+    the same "{+-}SCORE/DEPTH PV" convention match_play.py's
+    _format_eval_comment uses for self-play games (and
+    extract_selfplay_positions.py's COMMENT_RE expects on read) --
+    translates the engine's own raw "+MATE<N>"/"-MATE<N>" tellothers
+    token into that convention's "+M"/"-M" (no distance digit, matching
+    match_play.py exactly -- mate distance isn't independently verified
+    here either). The PV's own sentinel token (<EGTB>/<REC>/<=>/<#>),
+    when present, is passed through verbatim since that is exactly what
+    both extract_positions.py and extract_selfplay_positions.py key
+    their exactness detection off of."""
+    m = MATE_RE.match(score_tok)
+    if FLOAT_RE.match(score_tok):
+        score_str = score_tok
+    elif m:
+        score_str = "+M" if m.group(1) == "+" else "-M"
+    else:
+        score_str = score_tok  # unrecognized token: pass through as-is
+    comment = f"{score_str}/{depth}"
+    if pv:
+        comment += f" {pv}"
+    return comment
+
+
 def annotate_shard(shard_games_list, index, results):
     """Attach comments to each game's mainline nodes in place. index and
     results are both in script order and the same length; index maps
@@ -148,8 +188,8 @@ def annotate_shard(shard_games_list, index, results):
         while node.variations:
             node = node.variations[0]
             if pi in scores:
-                depth, score, pv = scores[pi]
-                node.comment = f"{score:+.2f}/{depth} {pv}"
+                depth, score_tok, pv = scores[pi]
+                node.comment = _format_comment(depth, score_tok, pv)
             pi += 1
 
 
diff --git a/annotate_pgns/build_training_tensors.py b/annotate_pgns/build_training_tensors.py
index 60e51e6..0f69d6c 100644
--- a/annotate_pgns/build_training_tensors.py
+++ b/annotate_pgns/build_training_tensors.py
@@ -12,26 +12,81 @@ index array + a row-offset array per perspective), exactly the input
 shape torch.nn.EmbeddingBag expects -- no fixed-width padding needed
 since bag sizes vary (3-32 active features per position).
 
-Skips exact rows (no p_search/r_outcome label -- see NNUE.md: a
-mate/EGTB/draw-claim score is a game-theoretic fact, not something to
-fit a sigmoid-WDL target against).
+Exact rows (mate/EGTB/draw-claim positions) are never fit against the
+sigmoid-WDL p_search target -- see NNUE.md: that score is a
+game-theoretic fact, not noisy evidence to calibrate a sigmoid
+against. But as of the EGTB-ground-truth-mixing decision, an exact row
+that extraction could assign a direct r_outcome (1.0/0.0/0.5,
+mover-relative) to is NOT dropped -- it is written to the tensors with
+target = r_outcome directly, lambda-blend bypassed entirely, which is
+exactly what the "mix TB ground truth into training" literature
+describes. Only rows with neither p_search nor a determinable
+r_outcome (unparseable comments, the "only"-move placeholder, an
+exact row from data predating this fix) are still skipped.
 
 Usage:
     python3 build_training_tensors.py \
         --in runs/sd10_twic_2M/games.train.labeled.quiet.jsonl \
         --outdir runs/sd10_twic_2M/tensors_train \
         --lambda 1.0
+
+Per-position feature-index computation (chess.Board(fen) + SEE-free
+active_feature_indices) is the expensive, embarrassingly-parallel part
+-- positions have no cross-dependency. --workers (default 16) shards
+the input across that many processes via multiprocessing.Pool,
+stitching results back together in input order so output is
+byte-for-byte identical to the old single-threaded path. --workers 1
+restores that original sequential path; --limit (smoke testing) always
+forces --workers 1 since it depends on exact file-order position
+counting.
 """
 import argparse
 import array
 import json
 import sys
+from multiprocessing import Pool
 from pathlib import Path
 
 import numpy as np
 
 from nnue_feature_spec import active_feature_indices, chess
 
+# Standard piece values for the auxiliary material-probe target
+# (nnue_model.py's material_probe, added 2026-10-17 after a measured
+# failure -- see that module's docstring). Plain material count,
+# mover-relative, NOT run through the WDL sigmoid -- this is a
+# straight regression target, deliberately non-saturating. Kings
+# excluded (always 1-1 on a legal board, carries no signal).
+_PIECE_VALUES = {
+    chess.PAWN: 1, chess.KNIGHT: 3, chess.BISHOP: 3,
+    chess.ROOK: 5, chess.QUEEN: 9,
+}
+
+
+def position_target(pos, lam):
+    """Returns the training target for one position, or None if it's
+    unusable (no p_search and no determinable r_outcome). Exact rows
+    with a direct r_outcome bypass p_search/lambda entirely -- see
+    this module's docstring and NNUE.md's EGTB-ground-truth-mixing
+    decision."""
+    if pos.get("exact"):
+        return pos.get("r_outcome")
+    if "p_search" not in pos:
+        return None
+    p_search = pos["p_search"]
+    r_outcome = pos.get("r_outcome", p_search)
+    return lam * p_search + (1.0 - lam) * r_outcome
+
+
+def material_balance(board):
+    """Mover-relative material balance in pawns, +ve = side to move
+    is materially ahead."""
+    total = 0
+    for piece_type, value in _PIECE_VALUES.items():
+        total += value * len(board.pieces(piece_type, chess.WHITE))
+        total -= value * len(board.pieces(piece_type, chess.BLACK))
+    return total if board.turn == chess.WHITE else -total
+
 
 def iter_positions(path):
     with open(path) as f:
@@ -41,20 +96,100 @@ def iter_positions(path):
                 yield pos
 
 
+def iter_chunks(iterable, chunk_size):
+    """Batches an iterable into lists of up to chunk_size items -- the
+    unit of work handed to one worker process at a time."""
+    chunk = []
+    for item in iterable:
+        chunk.append(item)
+        if len(chunk) >= chunk_size:
+            yield chunk
+            chunk = []
+    if chunk:
+        yield chunk
+
+
+def process_chunk(args):
+    """Runs in a worker process: the same per-position feature-index /
+    target / material computation main() used to do inline, just over
+    one chunk. Returns flat, already-concatenated arrays for the chunk
+    (not yet offset-adjusted -- that happens once, cheaply, back in
+    the main process when chunks are stitched together in order) plus
+    a per-position length list so the caller can rebuild the CSR
+    offset array without re-deriving it from the index array."""
+    positions, lam = args
+    white_idx = array.array("i")
+    black_idx = array.array("i")
+    white_lens = array.array("i")
+    black_lens = array.array("i")
+    stm = []
+    target = []
+    material = []
+    n_total = len(positions)
+    n_skipped = 0
+
+    for pos in positions:
+        t = position_target(pos, lam)
+        if t is None:
+            n_skipped += 1
+            continue
+
+        board = chess.Board(pos["fen"])
+        w = active_feature_indices(board, chess.WHITE)
+        b = active_feature_indices(board, chess.BLACK)
+
+        white_idx.extend(w)
+        white_lens.append(len(w))
+        black_idx.extend(b)
+        black_lens.append(len(b))
+
+        stm.append(1 if board.turn == chess.WHITE else 0)
+        target.append(t)
+        material.append(material_balance(board))
+
+    return (white_idx, white_lens, black_idx, black_lens, stm, target,
+            material, n_total, n_skipped)
+
+
 def main():
     ap = argparse.ArgumentParser(description=__doc__)
     ap.add_argument("--in", dest="infile", required=True)
     ap.add_argument("--outdir", required=True)
-    ap.add_argument("--lambda", dest="lam", type=float, default=1.0,
+    ap.add_argument("--lambda", dest="lam", type=float, default=0.9,
                      help="target = lambda*p_search + (1-lambda)*r_outcome "
-                          "(Phase A default 1.0, pure p_search -- see NNUE.md)")
+                          "for non-exact rows (exact rows always use "
+                          "r_outcome directly, ignoring lambda). Default "
+                          "0.9: with a Stockfish-sourced self-play corpus, "
+                          "p_search (SF's own WDL-calibrated search score) "
+                          "is a much lower-variance per-position signal "
+                          "than one game's actual outcome, so most of the "
+                          "weight stays on it -- see NNUE.md)")
     ap.add_argument("--limit", type=int, default=None,
-                     help="stop after this many usable positions (smoke testing)")
+                     help="stop after this many usable positions (smoke "
+                          "testing) -- forces --workers 1, since the limit "
+                          "is applied to the global used-position count in "
+                          "file order, which only a single sequential pass "
+                          "can honor exactly")
+    ap.add_argument("--workers", type=int, default=16,
+                     help="worker processes for the per-position feature-"
+                          "index computation (the expensive, embarrassingly"
+                          "-parallel part -- each position is independent). "
+                          "1 disables multiprocessing and runs the original "
+                          "sequential path")
+    ap.add_argument("--chunk-size", type=int, default=20000,
+                     help="positions per unit of work handed to a worker; "
+                          "only relevant when --workers > 1")
     args = ap.parse_args()
 
     outdir = Path(args.outdir)
     outdir.mkdir(parents=True, exist_ok=True)
 
+    workers = args.workers
+    if args.limit is not None and workers != 1:
+        print("--limit given: forcing --workers 1 for exact, "
+              "file-order-respecting smoke-test behavior", file=sys.stderr)
+        workers = 1
+
     # array.array('i') is a packed C int buffer -- avoids the ~28
     # bytes/element overhead of a Python list of boxed ints, which
     # matters once nnz reaches tens of millions.
@@ -64,39 +199,73 @@ def main():
     black_off = [0]
     stm = []       # 1 = white to move, 0 = black
     target = []
+    material = []  # auxiliary material-probe target, mover-relative pawns
 
     n_total = 0
     n_used = 0
     n_skipped = 0
 
-    for pos in iter_positions(args.infile):
-        n_total += 1
-        if args.limit is not None and n_used >= args.limit:
-            break
-        if pos.get("exact") or "p_search" not in pos:
-            n_skipped += 1
-            continue
+    if workers == 1:
+        for pos in iter_positions(args.infile):
+            n_total += 1
+            if args.limit is not None and n_used >= args.limit:
+                break
+            t = position_target(pos, args.lam)
+            if t is None:
+                n_skipped += 1
+                continue
 
-        board = chess.Board(pos["fen"])
-        w = active_feature_indices(board, chess.WHITE)
-        b = active_feature_indices(board, chess.BLACK)
+            board = chess.Board(pos["fen"])
+            w = active_feature_indices(board, chess.WHITE)
+            b = active_feature_indices(board, chess.BLACK)
 
-        white_idx.extend(w)
-        white_off.append(white_off[-1] + len(w))
-        black_idx.extend(b)
-        black_off.append(black_off[-1] + len(b))
+            white_idx.extend(w)
+            white_off.append(white_off[-1] + len(w))
+            black_idx.extend(b)
+            black_off.append(black_off[-1] + len(b))
 
-        stm.append(1 if board.turn == chess.WHITE else 0)
-        p_search = pos["p_search"]
-        r_outcome = pos.get("r_outcome", p_search)
-        target.append(args.lam * p_search + (1.0 - args.lam) * r_outcome)
+            stm.append(1 if board.turn == chess.WHITE else 0)
+            target.append(t)
+            material.append(material_balance(board))
+
+            n_used += 1
+            if n_used % 200000 == 0:
+                print(f"  ...{n_used} positions processed", file=sys.stderr)
+    else:
+        # Chunked across worker processes, but stitched back together
+        # in input order (Pool.imap preserves task order, unlike
+        # imap_unordered) -- output is byte-for-byte identical to the
+        # sequential path, just computed faster on an idle multi-core
+        # box. Each chunk's own CSR lengths get turned into this
+        # process's running cumulative offsets here, cheaply, since
+        # that's the one step that has to stay sequential.
+        chunks = ((chunk, args.lam)
+                  for chunk in iter_chunks(iter_positions(args.infile),
+                                            args.chunk_size))
+        with Pool(processes=workers) as pool:
+            for (w_idx, w_lens, b_idx, b_lens, chunk_stm, chunk_target,
+                 chunk_material, chunk_total,
+                 chunk_skipped) in pool.imap(process_chunk, chunks):
+                n_total += chunk_total
+                n_skipped += chunk_skipped
+
+                white_idx.extend(w_idx)
+                for length in w_lens:
+                    white_off.append(white_off[-1] + length)
+                black_idx.extend(b_idx)
+                for length in b_lens:
+                    black_off.append(black_off[-1] + length)
+
+                stm.extend(chunk_stm)
+                target.extend(chunk_target)
+                material.extend(chunk_material)
 
-        n_used += 1
-        if n_used % 200000 == 0:
-            print(f"  ...{n_used} positions processed", file=sys.stderr)
+                n_used += len(chunk_stm)
+                print(f"  ...{n_used} positions processed "
+                      f"({n_total} read)", file=sys.stderr)
 
     print(f"Read {n_total} rows ({args.infile})")
-    print(f"  {n_skipped} skipped (exact / no p_search label)")
+    print(f"  {n_skipped} skipped (no p_search and no determinable r_outcome)")
     print(f"  {n_used} usable positions")
 
     np.save(outdir / "white_idx.npy", np.frombuffer(white_idx, dtype=np.int32))
@@ -105,6 +274,7 @@ def main():
     np.save(outdir / "black_off.npy", np.array(black_off, dtype=np.int64))
     np.save(outdir / "stm.npy", np.array(stm, dtype=np.uint8))
     np.save(outdir / "target.npy", np.array(target, dtype=np.float32))
+    np.save(outdir / "material.npy", np.array(material, dtype=np.float32))
 
     nnz_white = len(white_idx)
     nnz_black = len(black_idx)
diff --git a/annotate_pgns/compare_vs_stockfish_net.py b/annotate_pgns/compare_vs_stockfish_net.py
new file mode 100644
index 0000000..0fe7514
--- /dev/null
+++ b/annotate_pgns/compare_vs_stockfish_net.py
@@ -0,0 +1,183 @@
+#!/usr/bin/env python3
+"""Compares our trained checkpoints against a real, matching-commit
+Stockfish build's NNUE evaluation on sampled real positions.
+
+NOT a weight-reuse step (see conversation: the goal is a from-scratch
+net) -- this is purely a diagnostic to see how far our net's score
+scale/ranking diverges from a known-strong reference, and whether
+feature factorization measurably closes that gap vs. the unfactored
+run.
+
+Stockfish's "NNUE evaluation" line (not "Final evaluation", which
+blends in a classical/hybrid term we have no equivalent of) is the
+directly comparable quantity -- both it and our forward() are pure
+NNUE output, no classical blending, both white-side-relative after
+sign conversion.
+
+Usage:
+    python3 compare_vs_stockfish_net.py \
+        --sample-file runs/stockfish_selfplay/games.holdout.labeled.quiet.jsonl \
+        --n-positions 1000 \
+        --sf-binary /tmp/sf_38e616/src/stockfish \
+        --sf-net nn-fdc1d0fe6455.nnue \
+        --ckpt-unfactored runs/stockfish_selfplay/nnue_ckpt_lambda075/ckpt_epoch033.pt \
+        --ckpt-factored runs/stockfish_selfplay/nnue_ckpt_lambda075_factored/ckpt_epoch033.pt
+"""
+import argparse
+import json
+import random
+import re
+import subprocess
+import sys
+
+import chess
+import numpy as np
+import torch
+
+from nnue_feature_spec import active_feature_indices
+from nnue_model import NnueNet
+
+NNUE_EVAL_RE = re.compile(r"NNUE evaluation\s+([+-]?\d+\.\d+) \(white side\)")
+
+
+def sample_fens(sample_file, n, seed):
+    rng = random.Random(seed)
+    reservoir = []
+    with open(sample_file) as f:
+        seen = 0
+        for line in f:
+            game = json.loads(line)
+            for pos in game["positions"]:
+                if pos.get("exact"):
+                    continue
+                seen += 1
+                if len(reservoir) < n:
+                    reservoir.append(pos["fen"])
+                else:
+                    j = rng.randint(0, seen - 1)
+                    if j < n:
+                        reservoir[j] = pos["fen"]
+    return reservoir
+
+
+def load_model(ckpt_path):
+    ckpt = torch.load(ckpt_path, map_location="cpu", weights_only=True)
+    model = NnueNet()
+    missing, unexpected = model.load_state_dict(ckpt["model"], strict=False)
+    bad_missing = [k for k in missing if k != "ft_p.weight"]
+    if bad_missing or unexpected:
+        sys.exit(f"{ckpt_path}: load_state_dict mismatch -- "
+                  f"missing={bad_missing}, unexpected={unexpected}")
+    model.eval()
+    return model
+
+
+def our_score_white_relative(model, fen):
+    board = chess.Board(fen)
+    w = active_feature_indices(board, chess.WHITE)
+    b = active_feature_indices(board, chess.BLACK)
+    w_idx = torch.tensor(w, dtype=torch.long)
+    w_off = torch.tensor([0], dtype=torch.long)
+    b_idx = torch.tensor(b, dtype=torch.long)
+    b_off = torch.tensor([0], dtype=torch.long)
+    stm = torch.tensor([1 if board.turn == chess.WHITE else 0], dtype=torch.long)
+    with torch.no_grad():
+        raw = model(w_idx, w_off, b_idx, b_off, stm).item()
+    return raw if board.turn == chess.WHITE else -raw
+
+
+def stockfish_scores_white_relative(sf_binary, sf_net, fens):
+    """One persistent Stockfish process, one 'position fen ...; eval'
+    round-trip per FEN, parsing the 'NNUE evaluation' line. Much
+    cheaper than spawning a fresh process per position."""
+    proc = subprocess.Popen(
+        [sf_binary], stdin=subprocess.PIPE, stdout=subprocess.PIPE,
+        stderr=subprocess.DEVNULL, text=True, bufsize=1)
+
+    def send(cmd):
+        proc.stdin.write(cmd + "\n")
+        proc.stdin.flush()
+
+    send(f"setoption name EvalFile value {sf_net}")
+    send("isready")
+    while True:
+        line = proc.stdout.readline()
+        if not line or "readyok" in line:
+            break
+
+    scores = []
+    for fen in fens:
+        send(f"position fen {fen}")
+        send("eval")
+        score = None
+        while True:
+            line = proc.stdout.readline()
+            if not line:
+                sys.exit("Stockfish process died mid-eval -- net likely "
+                         "failed to load (check EvalFile path/error output)")
+            m = NNUE_EVAL_RE.search(line)
+            if m:
+                score = float(m.group(1))
+            if "Final evaluation" in line:
+                break
+        scores.append(score)
+    send("quit")
+    proc.wait(timeout=5)
+    return scores
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--sample-file", required=True)
+    ap.add_argument("--n-positions", type=int, default=1000)
+    ap.add_argument("--seed", type=int, default=0)
+    ap.add_argument("--sf-binary", required=True)
+    ap.add_argument("--sf-net", required=True)
+    ap.add_argument("--ckpt-unfactored", required=True)
+    ap.add_argument("--ckpt-factored", required=True)
+    args = ap.parse_args()
+
+    print(f"Sampling {args.n_positions} non-exact positions from {args.sample_file} ...")
+    fens = sample_fens(args.sample_file, args.n_positions, args.seed)
+    print(f"  got {len(fens)} FENs")
+
+    print("Loading checkpoints ...")
+    model_unfactored = load_model(args.ckpt_unfactored)
+    model_factored = load_model(args.ckpt_factored)
+
+    print(f"Running {args.sf_binary} (net={args.sf_net}) ...")
+    sf_scores = stockfish_scores_white_relative(args.sf_binary, args.sf_net, fens)
+
+    print("Running our checkpoints ...")
+    ours_unfactored = [our_score_white_relative(model_unfactored, fen) for fen in fens]
+    ours_factored = [our_score_white_relative(model_factored, fen) for fen in fens]
+
+    sf = np.array(sf_scores, dtype=np.float64)
+    uf = np.array(ours_unfactored, dtype=np.float64)
+    fa = np.array(ours_factored, dtype=np.float64)
+
+    def report(name, ours):
+        diff = ours - sf
+        corr = np.corrcoef(ours, sf)[0, 1]
+        # best-fit scale (least squares, no intercept) -- how far off is
+        # our net's overall SCALE from Stockfish's, independent of the
+        # per-position noise captured by correlation/residual.
+        scale = float(np.dot(ours, sf) / np.dot(ours, ours))
+        resid_after_scale = ours * scale - sf
+        print(f"\n--- {name} vs Stockfish (n={len(sf)}) ---")
+        print(f"  correlation:            {corr:.4f}")
+        print(f"  raw mean|diff|:         {np.abs(diff).mean():.3f} pawns")
+        print(f"  raw median diff:        {np.median(diff):.3f} pawns")
+        print(f"  best-fit scale (ours->sf): {scale:.3f}")
+        print(f"  mean|resid| after scale-fit: {np.abs(resid_after_scale).mean():.3f} pawns")
+
+    report("Unfactored", uf)
+    report("Factored", fa)
+
+    np.savez("/tmp/compare_vs_sf.npz", sf=sf, unfactored=uf, factored=fa, fens=np.array(fens))
+    print("\nSaved raw arrays to /tmp/compare_vs_sf.npz")
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/annotate_pgns/export_nnue_weights.py b/annotate_pgns/export_nnue_weights.py
new file mode 100644
index 0000000..4c8cb99
--- /dev/null
+++ b/annotate_pgns/export_nnue_weights.py
@@ -0,0 +1,326 @@
+#!/usr/bin/env python3
+"""Exports a trained nnue_model.py checkpoint into a flat binary file
+the C engine (src/nnueeval.c's NnueLoadWeights) can read directly --
+no parsing, just a fixed-order sequence of arrays behind a small
+header.
+
+File format (all little-endian), version "NNU3":
+    magic:     4 bytes, ASCII "NNU3"
+    dims:      4x uint32 -- NNUE_DIMENSIONS, L1, L2, L3 (sanity-checked
+               against the C side's own compile-time constants on
+               load, so a mismatched build can't silently misread the
+               file)
+    ft_scale:  1x float32 -- fixed-point scale for ft.weight/ft_bias
+               (int16). NnueForward divides the int32 accumulator sum
+               by this to get the float "fc1 input" value.
+    act_scale: 3x float32 (AS1, AS2, AS3) -- fixed-point scales used
+               to quantize each stage's post-ReLU activation to uint8
+               before feeding the next layer's int8 dot product:
+                 AS1 quantizes fc1's input (the ReLU'd, concatenated
+                     accumulator halves)
+                 AS2 quantizes fc1's ReLU'd output (fc2's input)
+                 AS3 quantizes fc2's ReLU'd output (fc3's input)
+               Derived from real sampled activations (see
+               CalibrateActivationScales) with headroom, NOT from a
+               theoretical max -- see docs/NNUE.md's "Quantization
+               plan" / the int8 float-vs-quantized correctness check.
+    w_scale:   3x float32 (WS1, WS2, WS3) -- fixed-point scales for
+               fc1.weight/fc2.weight/fc3.weight (int8), one per layer,
+               each derived from that layer's own max|.| (exact, not
+               sampled -- weights are fixed, known precisely from the
+               checkpoint).
+    then, in this exact order:
+        ft.weight   [NNUE_DIMENSIONS x L1]   int16  (quantized @ ft_scale)
+        ft_bias     [L1]                     int16  (quantized @ ft_scale)
+        fc1.weight  [L2 x (2*L1)]            int8   (quantized @ WS1)
+        fc1.bias    [L2]                     float32 (unquantized --
+                                              added post-dequant, tiny
+                                              array, not worth it)
+        fc2.weight  [L3 x L2]                int8   (quantized @ WS2)
+        fc2.bias    [L3]                     float32
+        fc3.weight  [1 x L3]                 int8   (quantized @ WS3)
+        fc3.bias    [1]                      float32
+
+Both the embedding table (int16) and the FC head (int8) are quantized
+now -- see docs/NNUE.md's "Quantization plan". Measured accuracy cost
+of the FC-head int8 step, 5000 real holdout positions vs the exact
+float reference: mean |delta| 1.55cp, median 0.99cp, max 20.91cp, zero
+clipping -- an order of magnitude below the sd10-vs-sd16 search-depth
+noise (mean 19cp) the project already tolerates in its training
+labels. See the conversation/commit history for the simulation that
+produced those numbers; CalibrateActivationScales below reproduces the
+scale derivation deterministically from a real checkpoint + its
+training tensors.
+
+Usage:
+    python3 export_nnue_weights.py \
+        runs/sd10_twic_2M/nnue_ckpt/ckpt_epoch019.pt \
+        runs/sd10_twic_2M/nnue_ckpt/weights.nnue \
+        --tensor-dir runs/sd10_twic_2M/tensors_train
+"""
+import argparse
+import struct
+import sys
+from pathlib import Path
+
+import numpy as np
+import torch
+
+from nnue_feature_spec import DIMENSIONS, PS_NB
+from nnue_model import L1, L2, L3, NnueNet
+
+# Fixed-point scale for ft.weight/ft_bias. Checked against
+# ckpt_epoch019.pt (ft.weight max|.|=5.01, ft_bias max|.|=0.52): at
+# this scale, 5.01*4096 ~= 20502, comfortably under int16's 32767
+# ceiling (~1.6x headroom for a larger Phase B net before this needs
+# re-deriving). See docs/NNUE.md's "Quantization plan".
+FT_SCALE = 4096.0
+
+# Headroom multiplier applied to the largest value seen in a real
+# sample before deriving each activation scale -- protects against
+# rare tail values in the full corpus that a finite sample might not
+# have hit. 1.3x chosen as a round, conservative margin; revisit if a
+# larger/differently-distributed Phase B net ever reports clipping.
+ACT_HEADROOM = 1.3
+
+
+def CalibrateActivationScales(model, tensor_dir, n_sample=50000, seed=42):
+    """Runs n_sample real positions from tensor_dir (one of
+    build_training_tensors.py's output directories) through the float
+    model, capturing each stage's real post-ReLU activation
+    distribution via forward hooks, and returns (AS1, AS2, AS3) --
+    uint8 fixed-point scales with ACT_HEADROOM margin over the
+    observed max. Must run against the SAME checkpoint being
+    exported -- these scales are not portable across checkpoints."""
+    d = Path(tensor_dir)
+    white_idx = np.load(d / "white_idx.npy"); white_off = np.load(d / "white_off.npy")
+    black_idx = np.load(d / "black_idx.npy"); black_off = np.load(d / "black_off.npy")
+    stm = np.load(d / "stm.npy")
+    n = len(stm)
+    rng = np.random.default_rng(seed)
+    rows = rng.choice(n, size=min(n_sample, n), replace=False)
+    rows.sort()
+
+    captured = {}
+    h1 = model.fc1.register_forward_hook(
+        lambda m, i, o: captured.setdefault("fc1_in", []).append(i[0].detach()))
+    h2 = model.fc2.register_forward_hook(
+        lambda m, i, o: captured.setdefault("fc2_in", []).append(i[0].detach()))
+    h3 = model.fc3.register_forward_hook(
+        lambda m, i, o: captured.setdefault("fc3_in", []).append(i[0].detach()))
+
+    with torch.no_grad():
+        for i in range(0, len(rows), 5000):
+            chunk = rows[i:i + 5000]
+            w_idx_list, w_off_list, b_idx_list, b_off_list = [], [], [], []
+            off_w, off_b = 0, 0
+            for r in chunk:
+                wi = white_idx[white_off[r]:white_off[r + 1]]
+                bi = black_idx[black_off[r]:black_off[r + 1]]
+                w_idx_list.append(wi); w_off_list.append(off_w); off_w += len(wi)
+                b_idx_list.append(bi); b_off_list.append(off_b); off_b += len(bi)
+            w_idx = torch.from_numpy(np.concatenate(w_idx_list).astype(np.int64))
+            w_off = torch.tensor(w_off_list, dtype=torch.int64)
+            b_idx = torch.from_numpy(np.concatenate(b_idx_list).astype(np.int64))
+            b_off = torch.tensor(b_off_list, dtype=torch.int64)
+            stm_t = torch.from_numpy(stm[chunk].astype(np.int64))
+            model(w_idx, w_off, b_idx, b_off, stm_t)
+
+    h1.remove(); h2.remove(); h3.remove()
+
+    scales = []
+    for key in ("fc1_in", "fc2_in", "fc3_in"):
+        t = torch.cat(captured[key]).float()
+        m = t.max().item()
+        scales.append(255.0 / (m * ACT_HEADROOM))
+    return tuple(scales)
+
+
+def QuantizeFt(arr):
+    """Rounds + clamps a float32 array to int16 at FT_SCALE. Returns
+    (int16 ndarray, number of values that clipped)."""
+    scaled = np.round(arr.astype(np.float64) * FT_SCALE)
+    clipped = int(np.count_nonzero((scaled < -32768) | (scaled > 32767)))
+    scaled = np.clip(scaled, -32768, 32767)
+    return scaled.astype(np.int16), clipped
+
+
+def QuantizeWeightS8(arr, scale):
+    """Rounds + clamps a float32 array to int8 at the given scale.
+    Returns (int8 ndarray, number of values that clipped)."""
+    scaled = np.round(arr.astype(np.float64) * scale)
+    clipped = int(np.count_nonzero((scaled < -128) | (scaled > 127)))
+    scaled = np.clip(scaled, -128, 127)
+    return scaled.astype(np.int8), clipped
+
+
+def CoalesceFactorization(sd, tensor_dir, n_check=4096, seed=123):
+    """If the checkpoint has a trained ft_p factorization table (see
+    nnue_model.py's 2026-10-18 note), folds it into ft.weight in
+    place (idx % PS_NB is the P-factor's row for real feature idx)
+    and returns a modified copy of sd with ft_p dropped -- the
+    coalesced ft alone must reproduce the factored forward pass, so
+    ft_p carries zero inference-time cost. No-op (returns sd
+    unchanged) if ft_p isn't present (an older, unfactored
+    checkpoint).
+
+    Required correctness check, per NNUE.md: a bad coalesce silently
+    produces a WORSE net (real rows learn less because the factor
+    absorbed shared signal) with no sign of it in the training loss
+    curve. Verified here by comparing the pre-coalesce (factored,
+    two-table) forward pass against the post-coalesce (plain,
+    ft-only) forward pass on n_check real positions sampled from
+    tensor_dir -- must agree to ~1e-5 or this refuses to export."""
+    if "ft_p.weight" not in sd:
+        return sd
+
+    ft_p_weight = sd["ft_p.weight"]
+    if tuple(ft_p_weight.shape) != (PS_NB, L1):
+        sys.exit(f"ft_p.weight: expected shape {(PS_NB, L1)}, checkpoint "
+                  f"has {tuple(ft_p_weight.shape)} -- wrong model version?")
+
+    idx_mod = torch.arange(DIMENSIONS) % PS_NB
+    coalesced_ft_weight = sd["ft.weight"] + ft_p_weight[idx_mod]
+
+    sd_factored = dict(sd)
+    model_factored = NnueNet()
+    model_factored.load_state_dict(sd_factored)
+    model_factored.eval()
+
+    sd_plain = dict(sd)
+    sd_plain["ft.weight"] = coalesced_ft_weight
+    sd_plain["ft_p.weight"] = torch.zeros_like(ft_p_weight)
+    model_plain = NnueNet()
+    model_plain.load_state_dict(sd_plain)
+    model_plain.eval()
+
+    d = Path(tensor_dir)
+    white_idx = np.load(d / "white_idx.npy"); white_off = np.load(d / "white_off.npy")
+    black_idx = np.load(d / "black_idx.npy"); black_off = np.load(d / "black_off.npy")
+    stm = np.load(d / "stm.npy")
+    rng = np.random.default_rng(seed)
+    rows = rng.choice(len(stm), size=min(n_check, len(stm)), replace=False)
+    rows.sort()
+
+    w_idx_list, w_off_list, b_idx_list, b_off_list = [], [], [], []
+    off_w, off_b = 0, 0
+    for r in rows:
+        wi = white_idx[white_off[r]:white_off[r + 1]]
+        bi = black_idx[black_off[r]:black_off[r + 1]]
+        w_idx_list.append(wi); w_off_list.append(off_w); off_w += len(wi)
+        b_idx_list.append(bi); b_off_list.append(off_b); off_b += len(bi)
+    w_idx = torch.from_numpy(np.concatenate(w_idx_list).astype(np.int64))
+    w_off = torch.tensor(w_off_list, dtype=torch.int64)
+    b_idx = torch.from_numpy(np.concatenate(b_idx_list).astype(np.int64))
+    b_off = torch.tensor(b_off_list, dtype=torch.int64)
+    stm_t = torch.from_numpy(stm[rows].astype(np.int64))
+
+    with torch.no_grad():
+        out_factored = model_factored(w_idx, w_off, b_idx, b_off, stm_t)
+        out_plain = model_plain(w_idx, w_off, b_idx, b_off, stm_t)
+    max_delta = (out_factored - out_plain).abs().max().item()
+    if max_delta > 1e-4:
+        sys.exit(f"Coalesce check FAILED: max |factored - coalesced| = "
+                  f"{max_delta:.6g} over {n_check} real positions (expected "
+                  f"~1e-5) -- refusing to export a bad coalesce. Check "
+                  f"ft_p's shape/init and the idx % PS_NB formula.")
+    print(f"Coalesce check OK: max |factored - coalesced| = {max_delta:.3g} "
+          f"over {n_check} real positions")
+
+    sd_out = dict(sd)
+    sd_out["ft.weight"] = coalesced_ft_weight
+    del sd_out["ft_p.weight"]
+    return sd_out
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("checkpoint")
+    ap.add_argument("outfile")
+    ap.add_argument("--tensor-dir", required=True,
+                     help="a build_training_tensors.py output dir, used "
+                          "to calibrate the FC head's activation scales "
+                          "against real positions (e.g. tensors_train)")
+    args = ap.parse_args()
+
+    ckpt = torch.load(args.checkpoint, map_location="cpu", weights_only=True)
+    sd = ckpt["model"]
+    sd = CoalesceFactorization(sd, args.tensor_dir)
+
+    expected_shapes = {
+        "ft.weight": (DIMENSIONS, L1),
+        "ft_bias": (L1,),
+        "fc1.weight": (L2, 2 * L1),
+        "fc1.bias": (L2,),
+        "fc2.weight": (L3, L2),
+        "fc2.bias": (L3,),
+        "fc3.weight": (1, L3),
+        "fc3.bias": (1,),
+    }
+    for name, shape in expected_shapes.items():
+        if tuple(sd[name].shape) != shape:
+            sys.exit(f"{name}: expected shape {shape}, checkpoint has "
+                      f"{tuple(sd[name].shape)} -- wrong model version?")
+
+    model = NnueNet()
+    # strict=False: CoalesceFactorization deliberately drops ft_p.weight
+    # (folded into ft.weight above) when the checkpoint had it, but
+    # NnueNet always declares the ft_p parameter regardless -- missing_keys
+    # is expected to be exactly ["ft_p.weight"] in that case, or empty for
+    # an older unfactored checkpoint CoalesceFactorization no-op'd on.
+    missing, unexpected = model.load_state_dict(sd, strict=False)
+    bad_missing = [k for k in missing if k != "ft_p.weight"]
+    if bad_missing or unexpected:
+        sys.exit(f"load_state_dict mismatch beyond the expected ft_p.weight "
+                  f"drop -- missing={bad_missing}, unexpected={unexpected}")
+    model.eval()
+
+    AS1, AS2, AS3 = CalibrateActivationScales(model, args.tensor_dir)
+
+    ft_weight_f32 = sd["ft.weight"].detach().contiguous().to(torch.float32).numpy()
+    ft_bias_f32 = sd["ft_bias"].detach().contiguous().to(torch.float32).numpy()
+    fc1_w_f32 = sd["fc1.weight"].detach().contiguous().to(torch.float32).numpy()
+    fc2_w_f32 = sd["fc2.weight"].detach().contiguous().to(torch.float32).numpy()
+    fc3_w_f32 = sd["fc3.weight"].detach().contiguous().to(torch.float32).numpy()
+
+    ft_weight_q, clipped_ftw = QuantizeFt(ft_weight_f32)
+    ft_bias_q, clipped_ftb = QuantizeFt(ft_bias_f32)
+
+    WS1 = 127.0 / np.abs(fc1_w_f32).max()
+    WS2 = 127.0 / np.abs(fc2_w_f32).max()
+    WS3 = 127.0 / np.abs(fc3_w_f32).max()
+
+    fc1_w_q, clipped_w1 = QuantizeWeightS8(fc1_w_f32, WS1)
+    fc2_w_q, clipped_w2 = QuantizeWeightS8(fc2_w_f32, WS2)
+    fc3_w_q, clipped_w3 = QuantizeWeightS8(fc3_w_f32, WS3)
+
+    total_clipped = (clipped_ftw + clipped_ftb + clipped_w1 + clipped_w2 +
+                      clipped_w3)
+    if total_clipped:
+        sys.exit(f"Quantize: {total_clipped} value(s) clipped -- "
+                  f"re-check headroom/scales before shipping this export "
+                  f"(see docs/NNUE.md)")
+
+    with open(args.outfile, "wb") as f:
+        f.write(b"NNU3")
+        f.write(struct.pack("<4I", DIMENSIONS, L1, L2, L3))
+        f.write(struct.pack("<f", FT_SCALE))
+        f.write(struct.pack("<3f", AS1, AS2, AS3))
+        f.write(struct.pack("<3f", WS1, WS2, WS3))
+        f.write(ft_weight_q.tobytes())
+        f.write(ft_bias_q.tobytes())
+        f.write(fc1_w_q.tobytes())
+        f.write(sd["fc1.bias"].detach().contiguous().to(torch.float32).numpy().tobytes())
+        f.write(fc2_w_q.tobytes())
+        f.write(sd["fc2.bias"].detach().contiguous().to(torch.float32).numpy().tobytes())
+        f.write(fc3_w_q.tobytes())
+        f.write(sd["fc3.bias"].detach().contiguous().to(torch.float32).numpy().tobytes())
+
+    print(f"Wrote {args.outfile} from {args.checkpoint}")
+    print(f"  ft_scale={FT_SCALE}")
+    print(f"  act_scales (AS1,AS2,AS3)=({AS1:.4f},{AS2:.4f},{AS3:.4f})")
+    print(f"  weight_scales (WS1,WS2,WS3)=({WS1:.4f},{WS2:.4f},{WS3:.4f})")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/annotate_pgns/export_nnue_weights_v1_baseline.py b/annotate_pgns/export_nnue_weights_v1_baseline.py
new file mode 100644
index 0000000..57ac7e7
--- /dev/null
+++ b/annotate_pgns/export_nnue_weights_v1_baseline.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+"""One-off: re-create the OLD plain-float32 "NNU1" export format for a
+same-moment, apples-to-apples bench comparison against the new int16
+"NNU2" quantized format. Not the real exporter (export_nnue_weights.py
+is the real one now) -- throwaway, delete after use.
+"""
+import struct
+import sys
+import torch
+from nnue_feature_spec import DIMENSIONS
+from nnue_model import L1, L2, L3
+
+ckpt = torch.load(sys.argv[1], map_location="cpu", weights_only=True)
+sd = ckpt["model"]
+names = ["ft.weight", "ft_bias", "fc1.weight", "fc1.bias",
+         "fc2.weight", "fc2.bias", "fc3.weight", "fc3.bias"]
+with open(sys.argv[2], "wb") as f:
+    f.write(b"NNU1")
+    f.write(struct.pack("<4I", DIMENSIONS, L1, L2, L3))
+    for name in names:
+        arr = sd[name].detach().contiguous().to(torch.float32).numpy()
+        f.write(arr.tobytes())
+print(f"wrote {sys.argv[2]}")
diff --git a/annotate_pgns/extract_positions.py b/annotate_pgns/extract_positions.py
index 539419a..507f956 100644
--- a/annotate_pgns/extract_positions.py
+++ b/annotate_pgns/extract_positions.py
@@ -132,6 +132,19 @@ def main():
                 }
                 if reason:
                     rec["exact_reason"] = reason
+                # Mover-relative exact outcome, straight from the sign
+                # of the already-mover-relative score (unlike
+                # extract_selfplay_positions.py's foreign-engine path,
+                # typhoon's own "score" field here is never nulled out
+                # for mate/EGTB/REC rows -- only "only" (unsearched) and
+                # genuinely unparseable tokens carry no usable score).
+                # draw_recognizer ("<=>") is 0.5 regardless of sign.
+                if exact and reason == "draw_recognizer":
+                    rec["r_outcome"] = 0.5
+                elif exact and score is not None and reason not in (
+                        "unsearched_only_move",) and not (
+                        reason and reason.startswith("unrecognized")):
+                    rec["r_outcome"] = 1.0 if score > 0 else 0.0
                 out.write(json.dumps(rec) + "\n")
                 total += 1
                 if exact:
diff --git a/annotate_pgns/extract_selfplay_positions.py b/annotate_pgns/extract_selfplay_positions.py
new file mode 100644
index 0000000..5682b41
--- /dev/null
+++ b/annotate_pgns/extract_selfplay_positions.py
@@ -0,0 +1,247 @@
+#!/usr/bin/env python3
+"""Extract game-shaped records (game_id, result, positions[]) directly
+from a self-play PGN written by eval_tune/match_play.py's
+build_pgn_game -- a different source format from extract_positions.py
+(which mines annotate.py's raw shard_NN.script/.out logs), so this is
+a sibling script, not a variant of that one: the per-move score here
+already lives as a PGN comment on each move
+("{+-}NN.NN/DEPTH PV..." or "{+-}M/DEPTH PV..." for mate), not a
+tellothers log line, and needs a PGN parser (python-chess) rather than
+a log-line regex. Output schema matches build_games.py's games.jsonl
+exactly (that script's own docstring anticipated this: "a future
+loader for self-play games can produce the same games.jsonl schema
+and both sources mix freely at training-set-assembly time") -- drops
+straight into fit_k_and_label.py -> quiet_filter_games.py ->
+build_training_tensors.py unchanged, no separate split_holdout.py/
+build_games.py step needed since python-chess already gives real
+per-game grouping directly.
+
+2026-10-16, see docs/NNUE.md and the same day's qsearch-futility/NNUE-
+calibration investigation: this is the tool for pulling trainable
+positions out of the self-play corpus (pgn/self_play_games.pgn)
+specifically to fill the "decisive material, not yet mate" gap the
+TWIC/human-game corpus structurally can't supply (strong humans
+resign long before that point; these self-play games, generated with
+ResignThreshold disabled -- see match_play.py -- are deliberately
+played to real checkmate/stalemate/50-move/EGTB conclusions instead).
+
+Score sign convention matches extract_positions.py's output: each
+row's "score" is relative to the side to move AT "fen" (not a fixed
+White-POV), matching the training pipeline's existing convention
+elsewhere. match_play.py's own comments are stored White-relative
+(_format_eval_comment's score_white) -- converted back to mover-
+relative here.
+
+Per-move "mover" field records exactly which PGN White/Black header
+name played that move (e.g. "candidate", "baseline",
+"candidate_nnue_loose") -- lets a caller filter out any evaluator
+whose own self-reported scores shouldn't be trusted as training
+labels (see relabel_selfplay_pgn.py and the --exclude-movers flag
+below) without needing to know which games/dates those came from.
+
+Usage:
+    python3 extract_selfplay_positions.py \\
+        ../pgn/self_play_games.pgn \\
+        --out runs/selfplay_1013_1014/games.jsonl \\
+        --dates 2026-10-13 2026-10-14 \\
+        --event "typhoon self-play match" \\
+        --exclude-movers candidate_nnue_loose
+"""
+import argparse
+import json
+import re
+import sys
+
+import chess
+import chess.pgn
+
+COMMENT_RE = re.compile(
+    r"^([+-]M|[+-]?\d+\.\d+)(?:/(\d+))?(?:\s+(.*))?$"
+)
+
+# match_play.py's own mate-token path ("+M"/"-M", _format_eval_comment)
+# only fires when |score_white| >= NMATE_THRESHOLD=32000 (centipawns),
+# a threshold keyed to typhoon's own NMATE sentinel. It is never reached
+# when the mover is a *foreign* engine driven through a UCI<->xboard
+# adapter (see src/stockfish_xboard -- Stockfish via PolyGlot), because
+# PolyGlot's own UCI mate-score-to-xboard-cp conversion caps out far
+# below that: empirically, mate-in-N on this corpus (35,295
+# Stockfish-vs-Stockfish games, pgn/stockfish_self_play_training_data.pgn)
+# always lands as a plain float in [+90.00, +99.99] (or the negated
+# range), one per decisive game at depth 1 -- 5528 "+99.99/1" rows vs
+# exactly 5528 "1-0" results, 2728 "-99.99/1" vs 2729 "0-1". A direct
+# histogram of |score| confirms a clean gap: every value is either
+# below 70.00 or in [99.00, 99.99] (51,488 rows total, matching the
+# count of |score| >= 90 exactly) -- no ambiguous values near the
+# boundary. So any foreign-engine corpus needs this magnitude check in
+# addition to the literal "+M"/"-M" token case; MATE_CLAMP_THRESHOLD is
+# set well below the observed 99.00 floor (and well above any real
+# positional/material evaluation) to catch this with margin.
+MATE_CLAMP_THRESHOLD = 90.0
+
+
+def parse_comment(comment):
+    """Returns (score_white_pawns_or_None, depth_or_None, pv_or_None,
+    exact, exact_reason, outcome_white_or_None). score_white is None
+    for mate ("+M"/"-M", distance not recoverable from this comment
+    format -- see match_play.py's _format_eval_comment) -- exact=True
+    covers that case the same way extract_positions.py treats MATE
+    tokens, just without a distance number to carry. Also None for a
+    foreign-engine mate score that was clamped into an ordinary-looking
+    float instead of using match_play.py's own "+M"/"-M" convention
+    (see MATE_CLAMP_THRESHOLD above) -- same treatment, different
+    detection path.
+
+    outcome_white is the exact, White-relative game-theoretic result
+    (1.0/0.0/0.5) for exact rows where it can be read straight off this
+    comment, independent of score_white (which is deliberately nulled
+    out above for the mate/clamped cases -- the sign is still known
+    even though the magnitude isn't, and that sign is exactly the
+    label this field exists to carry). None when the comment doesn't
+    determine an outcome on its own (ordinary non-exact rows, and the
+    "unparseable" case) -- see NNUE.md's EGTB-ground-truth-mixing
+    decision: this is deliberately NOT derived from the game's final
+    result, since an exact-at-this-position verdict (TB/mate-forced)
+    is a fact about the position, not about whether the game later
+    actually reached that forced outcome."""
+    m = COMMENT_RE.match(comment.strip())
+    if not m:
+        return None, None, None, True, f"unparseable:{comment!r}", None
+    score_tok, depth_tok, pv = m.groups()
+    depth = int(depth_tok) if depth_tok else None
+    pv = pv.strip() if pv else None
+
+    exact, reason = False, None
+    if pv:
+        last_move = pv.rsplit(" ", 1)[-1]
+        if last_move in ("<EGTB>", "<REC>", "<=>", "<#>"):
+            exact = True
+            reason = {"<=>": "draw_recognizer",
+                      "<#>": "mate_recognizer"}.get(
+                          last_move, last_move.strip("<>").lower())
+
+    if reason == "draw_recognizer":
+        return None, depth, pv, True, reason, 0.5
+
+    if score_tok in ("+M", "-M"):
+        outcome_white = 1.0 if score_tok == "+M" else 0.0
+        return None, depth, pv, True, (reason or "mate"), outcome_white
+
+    score_white = float(score_tok)
+    if abs(score_white) >= MATE_CLAMP_THRESHOLD:
+        outcome_white = 1.0 if score_white > 0 else 0.0
+        return None, depth, pv, True, (reason or "mate_clamped"), outcome_white
+
+    if reason == "mate_recognizer":
+        outcome_white = 1.0 if score_white > 0 else 0.0
+        return score_white, depth, pv, True, reason, outcome_white
+
+    return score_white, depth, pv, exact, reason, None
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("pgn_file")
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--event", default=None,
+                     help="only extract games with this Event header "
+                          "(default: any)")
+    ap.add_argument("--dates", nargs="*", default=None,
+                     help="only extract games with Date header in this "
+                          "list (default: any)")
+    ap.add_argument("--exclude-movers", nargs="*", default=[],
+                     help="skip any ply whose mover's PGN name is in "
+                          "this list (e.g. candidate_nnue_loose) -- the "
+                          "move/FEN sequence for the rest of that same "
+                          "game is still usable, only that one side's "
+                          "own score annotations are dropped")
+    args = ap.parse_args()
+
+    target_dates = set(args.dates) if args.dates else None
+    exclude_movers = set(args.exclude_movers)
+
+    total_games = 0
+    games_in_scope = 0
+    total_rows = 0
+    total_exact = 0
+    total_excluded_mover = 0
+
+    with open(args.pgn_file, encoding="utf-8", errors="surrogateescape") as fin, \
+         open(args.out, "w") as fout:
+        while True:
+            game = chess.pgn.read_game(fin)
+            if game is None:
+                break
+            total_games += 1
+            if args.event is not None and game.headers.get("Event") != args.event:
+                continue
+            if target_dates is not None and game.headers.get("Date") not in target_dates:
+                continue
+            games_in_scope += 1
+
+            white_name = game.headers.get("White", "")
+            black_name = game.headers.get("Black", "")
+            game_result = game.headers.get("Result", "*")
+            board = game.board()
+            node = game
+            game_positions = []
+            ply_num = 0
+            while node.variations:
+                next_node = node.variations[0]
+                move = next_node.move
+                mover_is_white = board.turn == chess.WHITE
+                mover_name = white_name if mover_is_white else black_name
+                fen_before = board.fen()
+                comment = next_node.comment
+                board.push(move)
+                node = next_node
+                ply_num += 1
+                if not comment:
+                    continue
+                if mover_name in exclude_movers:
+                    total_excluded_mover += 1
+                    continue
+                (score_white, depth, pv, exact, reason,
+                 outcome_white) = parse_comment(comment)
+                score = None
+                if score_white is not None:
+                    score = score_white if mover_is_white else -score_white
+                rec = {
+                    "fen": fen_before,
+                    "depth": depth,
+                    "score": score,
+                    "score_raw": comment,
+                    "pv": pv,
+                    "exact": exact,
+                    "mover": mover_name,
+                }
+                if reason:
+                    rec["exact_reason"] = reason
+                if outcome_white is not None:
+                    rec["r_outcome"] = (outcome_white if mover_is_white
+                                         else 1.0 - outcome_white)
+                game_positions.append(rec)
+                total_rows += 1
+                if exact:
+                    total_exact += 1
+
+            if not game_positions:
+                continue
+            game_rec = {
+                "game_id": f"selfplay_{games_in_scope}_{game.headers.get('Date', '?')}_{ply_num}ply",
+                "result": game_result,
+                "positions": game_positions,
+            }
+            fout.write(json.dumps(game_rec) + "\n")
+
+    print(f"{total_games} games scanned, {games_in_scope} in scope "
+          f"(event={args.event!r}, dates={sorted(target_dates) if target_dates else None})",
+          file=sys.stderr)
+    print(f"{total_rows} positions written ({total_exact} exact), "
+          f"{total_excluded_mover} plies skipped for excluded mover",
+          file=sys.stderr)
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/annotate_pgns/fit_k_and_label.py b/annotate_pgns/fit_k_and_label.py
index 9fd09b0..11b2928 100644
--- a/annotate_pgns/fit_k_and_label.py
+++ b/annotate_pgns/fit_k_and_label.py
@@ -19,6 +19,16 @@ holdout data would leak information into what's supposed to be an
 untouched generalization check). The same fitted K is then applied
 (never refit) to label holdout positions too.
 
+Both the calibration-pair collection (--train) and the labeling pass
+(--train/--holdout) are sharded across --workers processes via
+multiprocessing.Pool, each worker vectorizing its own chunk of games
+with numpy -- this file used to run ~80 full Python-loop passes over
+every pair during the ternary search (2 mse() evals x 40 iterations),
+which was the slowest step in the post-annotation pipeline relative to
+how little actual compute it does. Per-game-chunk numpy vectorization
+plus Pool.imap (order-preserving, like build_training_tensors.py's
+existing pattern) turns that into array ops.
+
 Usage:
     python3 fit_k_and_label.py --train runs/sd10_twic_2M/games.train.jsonl \
         --holdout runs/sd10_twic_2M/games.holdout.jsonl \
@@ -27,7 +37,9 @@ Usage:
 """
 import argparse
 import json
-import math
+from multiprocessing import Pool
+
+import numpy as np
 
 
 def result_to_white_score(result):
@@ -35,77 +47,160 @@ def result_to_white_score(result):
     return {"1-0": 1.0, "0-1": 0.0, "1/2-1/2": 0.5}.get(result)
 
 
-def stm_relative_outcome(white_score, fen):
-    """Flip a White-perspective game result to the side-to-move's own
-    perspective at this FEN, matching the score field's own
-    side-to-move-relative convention."""
-    stm_is_white = fen.split()[1] == "w"
-    return white_score if stm_is_white else (1.0 - white_score)
-
-
-def collect_calibration_pairs(games_path):
-    """(score, r_outcome) for every non-exact position across all
-    games in this file."""
-    pairs = []
-    n_unknown_result = 0
-    with open(games_path) as f:
-        for line in f:
-            g = json.loads(line)
-            white_score = result_to_white_score(g["result"])
-            if white_score is None:
-                n_unknown_result += 1
-                continue
-            for pos in g["positions"]:
-                if pos.get("exact"):
-                    continue
-                r = stm_relative_outcome(white_score, pos["fen"])
-                pairs.append((pos["score"], r))
-    return pairs, n_unknown_result
+def sigmoid(scores, k):
+    return 1.0 / (1.0 + np.exp(-scores / k))
 
 
-def mse(pairs, k):
-    err = 0.0
-    for score, r in pairs:
-        p = 1.0 / (1.0 + math.exp(-score / k))
-        err += (r - p) ** 2
-    return err / max(len(pairs), 1)
+def mse(scores, outcomes, k):
+    p = sigmoid(scores, k)
+    return float(np.mean((outcomes - p) ** 2)) if len(scores) else 0.0
 
 
-def fit_k(pairs, lo=0.05, hi=20.0, iterations=40):
+def fit_k(scores, outcomes, lo=0.05, hi=20.0, iterations=40):
     for _ in range(iterations):
         m1 = lo + (hi - lo) / 3
         m2 = hi - (hi - lo) / 3
-        if mse(pairs, m1) < mse(pairs, m2):
+        if mse(scores, outcomes, m1) < mse(scores, outcomes, m2):
             hi = m2
         else:
             lo = m1
     return (lo + hi) / 2
 
 
-def label_file(in_path, out_path, k):
+def iter_chunks(path, chunk_size):
+    """Batches raw lines (one game per line) of an open jsonl file into
+    lists of up to chunk_size lines -- the unit of work handed to one
+    worker process at a time. Text, not parsed JSON, so the parse
+    itself happens in the worker, same division of labor as
+    build_training_tensors.py's iter_chunks."""
+    chunk = []
+    with open(path) as f:
+        for line in f:
+            chunk.append(line)
+            if len(chunk) >= chunk_size:
+                yield chunk
+                chunk = []
+    if chunk:
+        yield chunk
+
+
+def _non_exact_scores_and_outcomes(game):
+    """(scores, outcomes) numpy arrays for one game's non-exact
+    positions, or None if the game's result is unknown or it has no
+    non-exact positions."""
+    white_score = result_to_white_score(game["result"])
+    if white_score is None:
+        return None
+    non_exact = [p for p in game["positions"] if not p.get("exact")]
+    if not non_exact:
+        return None
+    scores = np.fromiter((p["score"] for p in non_exact), dtype=np.float64,
+                          count=len(non_exact))
+    is_white = np.fromiter(
+        (p["fen"].split()[1] == "w" for p in non_exact), dtype=bool,
+        count=len(non_exact))
+    outcomes = np.where(is_white, white_score, 1.0 - white_score)
+    return scores, outcomes
+
+
+def collect_chunk(lines):
+    """Runs in a worker process: parses a chunk of game lines and
+    returns (scores, outcomes, n_unknown_result) for every non-exact
+    position across those games."""
+    score_pieces = []
+    outcome_pieces = []
+    n_unknown_result = 0
+    for line in lines:
+        g = json.loads(line)
+        result = _non_exact_scores_and_outcomes(g)
+        if result is None:
+            if result_to_white_score(g["result"]) is None:
+                n_unknown_result += 1
+            continue
+        scores, outcomes = result
+        score_pieces.append(scores)
+        outcome_pieces.append(outcomes)
+    if not score_pieces:
+        return np.empty(0), np.empty(0), n_unknown_result
+    return (np.concatenate(score_pieces), np.concatenate(outcome_pieces),
+            n_unknown_result)
+
+
+def collect_calibration_pairs(games_path, workers, chunk_size):
+    """(scores, outcomes) numpy arrays for every non-exact position
+    across all games in this file, built by sharding games across
+    --workers processes."""
+    score_pieces = []
+    outcome_pieces = []
+    n_unknown_result = 0
+    chunks = iter_chunks(games_path, chunk_size)
+    with Pool(processes=workers) as pool:
+        for scores, outcomes, n_unk in pool.imap(collect_chunk, chunks):
+            if len(scores):
+                score_pieces.append(scores)
+                outcome_pieces.append(outcomes)
+            n_unknown_result += n_unk
+    if not score_pieces:
+        return np.empty(0), np.empty(0), n_unknown_result
+    return (np.concatenate(score_pieces), np.concatenate(outcome_pieces),
+            n_unknown_result)
+
+
+def label_chunk(args):
+    """Runs in a worker process: labels every non-exact position in a
+    chunk of game lines (vectorized per game with numpy) and returns
+    the re-serialized lines plus counts, in the same order they came
+    in."""
+    lines, k = args
+    out_lines = []
     n_games = 0
     n_positions = 0
     n_exact_skipped = 0
     n_unknown_result = 0
-    with open(in_path) as fin, open(out_path, "w") as fout:
-        for line in fin:
-            g = json.loads(line)
-            white_score = result_to_white_score(g["result"])
-            if white_score is None:
-                n_unknown_result += 1
-                fout.write(json.dumps(g) + "\n")
-                continue
-            for pos in g["positions"]:
-                n_positions += 1
-                if pos.get("exact"):
-                    n_exact_skipped += 1
-                    continue
-                r = stm_relative_outcome(white_score, pos["fen"])
-                p_search = 1.0 / (1.0 + math.exp(-pos["score"] / k))
-                pos["p_search"] = p_search
-                pos["r_outcome"] = r
-            fout.write(json.dumps(g) + "\n")
-            n_games += 1
+    for line in lines:
+        g = json.loads(line)
+        white_score = result_to_white_score(g["result"])
+        positions = g["positions"]
+        n_positions += len(positions)
+        if white_score is None:
+            n_unknown_result += 1
+            out_lines.append(json.dumps(g))
+            continue
+
+        idxs = [i for i, p in enumerate(positions) if not p.get("exact")]
+        n_exact_skipped += len(positions) - len(idxs)
+        if idxs:
+            scores = np.fromiter((positions[i]["score"] for i in idxs),
+                                  dtype=np.float64, count=len(idxs))
+            is_white = np.fromiter(
+                (positions[i]["fen"].split()[1] == "w" for i in idxs),
+                dtype=bool, count=len(idxs))
+            outcomes = np.where(is_white, white_score, 1.0 - white_score)
+            p_search = sigmoid(scores, k)
+            for j, i in enumerate(idxs):
+                positions[i]["p_search"] = float(p_search[j])
+                positions[i]["r_outcome"] = float(outcomes[j])
+
+        out_lines.append(json.dumps(g))
+        n_games += 1
+    return out_lines, n_games, n_positions, n_exact_skipped, n_unknown_result
+
+
+def label_file(in_path, out_path, k, workers, chunk_size):
+    n_games = 0
+    n_positions = 0
+    n_exact_skipped = 0
+    n_unknown_result = 0
+    chunks = ((chunk, k) for chunk in iter_chunks(in_path, chunk_size))
+    with Pool(processes=workers) as pool, open(out_path, "w") as fout:
+        for (out_lines, chunk_games, chunk_positions, chunk_exact,
+             chunk_unknown) in pool.imap(label_chunk, chunks):
+            for line in out_lines:
+                fout.write(line + "\n")
+            n_games += chunk_games
+            n_positions += chunk_positions
+            n_exact_skipped += chunk_exact
+            n_unknown_result += chunk_unknown
     return n_games, n_positions, n_exact_skipped, n_unknown_result
 
 
@@ -115,26 +210,56 @@ def main():
     ap.add_argument("--holdout", required=True)
     ap.add_argument("--train-out", required=True)
     ap.add_argument("--holdout-out", required=True)
+    ap.add_argument("--workers", type=int, default=16,
+                     help="worker processes for both the calibration-pair "
+                          "collection and the labeling pass (default 16, "
+                          "kept below this box's 24 cores to leave other "
+                          "jobs some breathing room)")
+    ap.add_argument("--chunk-size", type=int, default=2000,
+                     help="games per unit of work handed to a worker")
+    ap.add_argument("--fixed-k", type=float, default=None,
+                     help="use this K directly instead of re-fitting it "
+                          "from --train. 2026-10-16: needed when labeling "
+                          "a small supplementary slice (e.g. a self-play "
+                          "augmentation batch) to continue training an "
+                          "EXISTING checkpoint -- that checkpoint's "
+                          "weights were fit against labels built with its "
+                          "own K (nnue_model.py's DEFAULT_K unless noted "
+                          "otherwise), and re-fitting K from a small, "
+                          "distributionally skewed slice (e.g. one "
+                          "deliberately overweighted toward decisive-"
+                          "material positions) would silently relabel "
+                          "everything on a different, inconsistent scale "
+                          "-- the exact mismatch this module's own closing "
+                          "warning describes. Omit for a normal from-"
+                          "scratch fit on a large, representative pool.")
     args = ap.parse_args()
 
-    print(f"Collecting calibration pairs from {args.train} (TRAIN only) ...")
-    pairs, n_unknown = collect_calibration_pairs(args.train)
-    print(f"  {len(pairs)} (score, outcome) pairs "
-          f"({n_unknown} games skipped: unknown result)")
+    if args.fixed_k is not None:
+        k = args.fixed_k
+        print(f"Using fixed K = {k:.4f} (not re-fit from {args.train!r})")
+    else:
+        print(f"Collecting calibration pairs from {args.train} (TRAIN only) ...")
+        scores, outcomes, n_unknown = collect_calibration_pairs(
+            args.train, args.workers, args.chunk_size)
+        print(f"  {len(scores)} (score, outcome) pairs "
+              f"({n_unknown} games skipped: unknown result)")
 
-    print("Fitting K (ternary search, minimizing sigmoid-MSE against "
-          "actual game outcomes) ...")
-    k = fit_k(pairs)
-    final_mse = mse(pairs, k)
-    print(f"  K = {k:.2f}  (MSE at fit = {final_mse:.5f})")
+        print("Fitting K (ternary search, minimizing sigmoid-MSE against "
+              "actual game outcomes) ...")
+        k = fit_k(scores, outcomes)
+        final_mse = mse(scores, outcomes, k)
+        print(f"  K = {k:.2f}  (MSE at fit = {final_mse:.5f})")
 
     print(f"\nLabeling {args.train} -> {args.train_out} ...")
-    n_g, n_p, n_exact, n_unk = label_file(args.train, args.train_out, k)
+    n_g, n_p, n_exact, n_unk = label_file(
+        args.train, args.train_out, k, args.workers, args.chunk_size)
     print(f"  {n_g} games, {n_p} positions, {n_exact} exact (unlabeled), "
           f"{n_unk} unknown-result games (unlabeled)")
 
     print(f"Labeling {args.holdout} -> {args.holdout_out} (same K, not refit) ...")
-    n_g, n_p, n_exact, n_unk = label_file(args.holdout, args.holdout_out, k)
+    n_g, n_p, n_exact, n_unk = label_file(
+        args.holdout, args.holdout_out, k, args.workers, args.chunk_size)
     print(f"  {n_g} games, {n_p} positions, {n_exact} exact (unlabeled), "
           f"{n_unk} unknown-result games (unlabeled)")
 
diff --git a/annotate_pgns/nnue_model.py b/annotate_pgns/nnue_model.py
index acdc8fd..c483ce7 100644
--- a/annotate_pgns/nnue_model.py
+++ b/annotate_pgns/nnue_model.py
@@ -16,11 +16,54 @@ activations / quantization-aware scaling (this is a float32-only
 Phase A prototype, no export-for-C-inference step exists yet), the
 lazy-eval PSQT-term fast path, and the super-lazy material-only
 prefilter. All irrelevant to proving the training pipeline works.
+
+2026-10-17: added `material_probe`, a small auxiliary linear head on
+the same concatenated accumulator `fc1` reads, trained jointly with
+the main WDL loss (see nnue_train.py's --material-loss-weight) against
+each position's raw material balance (mover-relative pawns, no
+sigmoid). Added after a real, measured failure: a net trained with
+plain MSE (and separately, exponent=2.6) against WDL-sigmoid targets
+alone could not tell a lost pawn from a lost queen in crowded midgame
+positions (see sanity_check_net.py) -- directionally correct but
+magnitude-blind. The probe supplies a non-saturating gradient straight
+at material, forcing `ft`'s embedding weights to keep material
+linearly recoverable even where the main WDL loss has little pressure
+to care. It's training-only scaffolding, like the SF doc's "virtual
+features" factorization trick -- dropped before export, zero
+inference-time or quantization-format cost.
+
+2026-10-18: added `ft_p`, the actual Stockfish-doc "P" factor (see
+NNUE.md's "feature factorization" status entry) -- a second,
+king-bucket-INDEPENDENT EmbeddingBag over PS_NB=704 rows
+(piece_square x piece_type x piece_color only), summed into the same
+accumulator `ft` feeds. Motivation: HalfKAv2_hm's real 22528-row
+table is king-bucket-conditioned, so a (piece, square, color) tuple
+that's common in general chess but rare in a specific king bucket
+gets its own mostly-untrained row -- measured directly on the
+Stockfish self-play training pool: 9.9% of rows never fire at all,
+23.6% see fewer than 50 examples. `ft_p`'s 704 rows see ~32x more
+observations each (one per king bucket collapsed into one), supplying
+dense gradient to a king-bucket-independent piece-value prior that
+`ft`'s sparser rows can then specialize away from. Zero-initialized
+(nn.init.zeros_) specifically so it contributes nothing at step 0 --
+a default (unit-variance) init would silently double the
+accumulator's effective scale at init, confounding the factorization
+effect with an incidental learning-rate change. Deliberately no bias
+of its own -- `ft_bias` stays the only bias on this sum, which is
+also what keeps the post-training coalescing step (see
+export_nnue_weights.py) exact rather than approximate. Training-only,
+like material_probe: coalesced into `ft.weight` and dropped before
+export, since `idx % PS_NB` (the P-factor's index for real feature
+`idx`) is recoverable from the real index alone, needing no change to
+build_training_tensors.py's cached tensors. Deliberately NOT adding
+Stockfish's other two documented factors ("K", "HalfRelativeKP") --
+HalfKAv2_hm already encodes the king directly via the shared PS_KING
+plane, making a K factor redundant here; one variable at a time.
 """
 import torch
 import torch.nn as nn
 
-from nnue_feature_spec import DIMENSIONS
+from nnue_feature_spec import DIMENSIONS, PS_NB
 
 L1 = 256   # feature-transformer output width (per perspective)
 L2 = 32
@@ -30,7 +73,8 @@ DEFAULT_K = 1.3123   # pawns; fit_k_and_label.py's Phase-A train-split fit
 
 
 class NnueNet(nn.Module):
-    def __init__(self, dimensions=DIMENSIONS, l1=L1, l2=L2, l3=L3):
+    def __init__(self, dimensions=DIMENSIONS, l1=L1, l2=L2, l3=L3,
+                 ps_nb=PS_NB):
         super().__init__()
         # Single shared feature transformer: one [dimensions, l1]
         # weight matrix, applied once per perspective to that
@@ -40,15 +84,37 @@ class NnueNet(nn.Module):
         self.ft = nn.EmbeddingBag(dimensions, l1, mode="sum")
         self.ft_bias = nn.Parameter(torch.zeros(l1))
 
+        # Training-only "P" factorization table -- see module
+        # docstring's 2026-10-18 note. King-bucket-independent
+        # (ps_nb=704 rows, vs ft's full 22528), summed into the same
+        # accumulator. Zero-init so it contributes nothing at step 0.
+        self.ft_p = nn.EmbeddingBag(ps_nb, l1, mode="sum")
+        nn.init.zeros_(self.ft_p.weight)
+        self.ps_nb = ps_nb
+
         self.fc1 = nn.Linear(2 * l1, l2)
         self.fc2 = nn.Linear(l2, l3)
         self.fc3 = nn.Linear(l3, 1)
 
-    def forward(self, white_idx, white_off, black_idx, black_off, stm):
-        """stm: float/long tensor of shape [B], 1.0 if white to move
-        else 0.0. Returns raw score in pawn-ish units, shape [B]."""
-        white_acc = self.ft(white_idx, white_off) + self.ft_bias
-        black_acc = self.ft(black_idx, black_off) + self.ft_bias
+        # Auxiliary, training-only head -- see module docstring's
+        # 2026-10-17 note. Reads the same concatenated accumulator
+        # fc1 does; not used at inference/export time.
+        self.material_probe = nn.Linear(2 * l1, 1)
+
+    def accumulator(self, white_idx, white_off, black_idx, black_off, stm):
+        """The ReLU'd, own-perspective-first concatenated accumulator
+        (shape [B, 2*l1]) -- shared input to both fc1 (the real eval
+        head) and material_probe (the auxiliary head). ft_p's
+        contribution uses the identical bag offsets as ft -- same
+        positions, same active-feature counts, just a second lookup
+        of each index modulo ps_nb into the smaller, denser table."""
+        white_idx_p = white_idx % self.ps_nb
+        black_idx_p = black_idx % self.ps_nb
+
+        white_acc = (self.ft(white_idx, white_off)
+                     + self.ft_p(white_idx_p, white_off) + self.ft_bias)
+        black_acc = (self.ft(black_idx, black_off)
+                     + self.ft_p(black_idx_p, black_off) + self.ft_bias)
 
         white_acc = torch.relu(white_acc)
         black_acc = torch.relu(black_acc)
@@ -57,13 +123,29 @@ class NnueNet(nn.Module):
         # own-perspective-first concatenation, per side to move.
         own = stm * white_acc + (1.0 - stm) * black_acc
         enemy = stm * black_acc + (1.0 - stm) * white_acc
-        x = torch.cat([own, enemy], dim=1)
+        return torch.cat([own, enemy], dim=1)
 
+    def forward(self, white_idx, white_off, black_idx, black_off, stm):
+        """stm: float/long tensor of shape [B], 1.0 if white to move
+        else 0.0. Returns raw score in pawn-ish units, shape [B]."""
+        x = self.accumulator(white_idx, white_off, black_idx, black_off, stm)
         x = torch.relu(self.fc1(x))
         x = torch.relu(self.fc2(x))
         x = self.fc3(x)
         return x.view(-1)
 
+    def forward_with_material(self, white_idx, white_off, black_idx, black_off, stm):
+        """Like forward(), but also returns the auxiliary material
+        probe's prediction (mover-relative pawns, no sigmoid) off the
+        same accumulator -- one extra Linear(2*l1, 1), training-only.
+        Returns (raw_score, material_pred), both shape [B]."""
+        x = self.accumulator(white_idx, white_off, black_idx, black_off, stm)
+        material_pred = self.material_probe(x).view(-1)
+        y = torch.relu(self.fc1(x))
+        y = torch.relu(self.fc2(y))
+        y = self.fc3(y)
+        return y.view(-1), material_pred
+
     def win_probability(self, white_idx, white_off, black_idx, black_off, stm, k=DEFAULT_K):
         raw = self.forward(white_idx, white_off, black_idx, black_off, stm)
         return torch.sigmoid(raw / k)
diff --git a/annotate_pgns/nnue_train.py b/annotate_pgns/nnue_train.py
index 6112ab3..ea3d0a9 100644
--- a/annotate_pgns/nnue_train.py
+++ b/annotate_pgns/nnue_train.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python3
 """Phase A training loop: local CPU, PyTorch, over the cached CSR
-tensors build_training_tensors.py produces. Goal is pipeline/
+stensors build_training_tensors.py produces. Goal is pipeline/
 framework validation (per NNUE.md) -- prove the data loader, model,
 loss, and checkpointing work end-to-end and get fast local iteration
 when something breaks, not to produce the best possible net. The
@@ -17,7 +17,7 @@ Usage:
         --train-dir runs/sd10_twic_2M/tensors_train \
         --holdout-dir runs/sd10_twic_2M/tensors_holdout \
         --outdir runs/sd10_twic_2M/nnue_ckpt \
-        --epochs 20 --batch-size 8192 --threads 8
+        --epochs 20 --batch-size 8192 --threads 16
 """
 import argparse
 import json
@@ -26,7 +26,6 @@ from pathlib import Path
 
 import numpy as np
 import torch
-import torch.nn as nn
 
 from nnue_model import NnueNet, DEFAULT_K
 
@@ -43,10 +42,17 @@ class CsrTensors:
         self.black_off = np.load(tensor_dir / "black_off.npy")
         self.stm = np.load(tensor_dir / "stm.npy")
         self.target = np.load(tensor_dir / "target.npy")
+        # material.npy added 2026-10-17 (see nnue_model.py's
+        # material_probe) -- optional for backward compat with
+        # tensor dirs built before this field existed.
+        material_path = tensor_dir / "material.npy"
+        self.material = (np.load(material_path) if material_path.exists()
+                          else np.zeros_like(self.target))
         self.n = len(self.target)
         assert len(self.white_off) == self.n + 1
         assert len(self.black_off) == self.n + 1
         assert len(self.stm) == self.n
+        assert len(self.material) == self.n
 
     def __len__(self):
         return self.n
@@ -68,6 +74,16 @@ def _gather_batch(idx_arr, off_arr, row_ids):
     return flat, batch_off
 
 
+def _to_device(batch, device):
+    """Moves one iterate_batches() tuple onto device. A no-op (same
+    tensors returned) when device is "cpu", since CsrTensors already
+    hands back CPU tensors directly -- only a CUDA run pays this
+    transfer cost."""
+    if device == "cpu":
+        return batch
+    return tuple(t.to(device, non_blocking=True) for t in batch)
+
+
 def iterate_batches(data: CsrTensors, batch_size, shuffle, rng):
     order = rng.permutation(data.n) if shuffle else np.arange(data.n)
     for start in range(0, data.n, batch_size):
@@ -77,6 +93,7 @@ def iterate_batches(data: CsrTensors, batch_size, shuffle, rng):
         b_idx, b_off = _gather_batch(data.black_idx, data.black_off, rows_sorted)
         stm = data.stm[rows_sorted]
         target = data.target[rows_sorted]
+        material = data.material[rows_sorted]
         yield (
             torch.from_numpy(w_idx.astype(np.int64, copy=False)),
             torch.from_numpy(w_off),
@@ -84,23 +101,51 @@ def iterate_batches(data: CsrTensors, batch_size, shuffle, rng):
             torch.from_numpy(b_off),
             torch.from_numpy(stm.astype(np.int64, copy=False)),
             torch.from_numpy(target.copy()),
+            torch.from_numpy(material.astype(np.float32, copy=True)),
         )
 
 
-def evaluate(model, data, batch_size, k):
+def powered_loss(pred, target, exponent, reduction="mean"):
+    """|pred - target|^exponent, exponent=2 is plain MSE. Stockfish's
+    own training notes (nnue-pytorch wiki) report good results with
+    2.6 -- widens the loss ratio between small and large WDL-space
+    errors versus plain MSE, pushing the optimizer harder to fix
+    decisively-wrong (e.g. material-blind) predictions relative to
+    subtle positional ones. 2026-10-17: tried after plain MSE (exp=2)
+    showed near-identical score drops whether a pawn or a queen was
+    removed from a position (see sanity_check_net.py)."""
+    err = (pred - target).abs().pow(exponent)
+    return err.sum() if reduction == "sum" else err.mean()
+
+
+def evaluate(model, data, batch_size, k, loss_exponent, device="cpu"):
+    """Returns (wdl_loss, material_mse) -- kept separate rather than
+    pre-combined so the log/checkpoint can show whether the material
+    probe is actually converging, independently of whatever weight
+    the caller gives it in the training loss."""
     model.eval()
-    total_loss = 0.0
+    total_wdl_loss = 0.0
+    total_material_loss = 0.0
     total_n = 0
     rng = np.random.default_rng(0)
     with torch.no_grad():
-        for w_idx, w_off, b_idx, b_off, stm, target in iterate_batches(
+        for w_idx, w_off, b_idx, b_off, stm, target, material in iterate_batches(
                 data, batch_size, shuffle=False, rng=rng):
-            pred = model.win_probability(w_idx, w_off, b_idx, b_off, stm, k=k)
-            loss = nn.functional.mse_loss(pred, target, reduction="sum")
-            total_loss += loss.item()
+            w_idx, w_off, b_idx, b_off, stm, target, material = _to_device(
+                (w_idx, w_off, b_idx, b_off, stm, target, material), device)
+            raw, material_pred = model.forward_with_material(w_idx, w_off, b_idx, b_off, stm)
+            pred = torch.sigmoid(raw / k)
+            wdl_loss = powered_loss(pred, target, loss_exponent, reduction="sum")
+            material_loss = nn_mse_sum(material_pred, material)
+            total_wdl_loss += wdl_loss.item()
+            total_material_loss += material_loss.item()
             total_n += len(target)
     model.train()
-    return total_loss / max(total_n, 1)
+    return total_wdl_loss / max(total_n, 1), total_material_loss / max(total_n, 1)
+
+
+def nn_mse_sum(pred, target):
+    return (pred - target).pow(2).sum()
 
 
 def main():
@@ -112,16 +157,65 @@ def main():
     ap.add_argument("--batch-size", type=int, default=8192)
     ap.add_argument("--lr", type=float, default=1e-3)
     ap.add_argument("--k", type=float, default=DEFAULT_K)
-    ap.add_argument("--threads", type=int, default=8,
-                     help="torch CPU thread count -- kept modest by default "
-                          "since this is a shared multi-user box")
+    ap.add_argument("--loss-exponent", type=float, default=2.0,
+                     help="exponent on |pred-target| in the loss (2.0 = "
+                          "plain MSE). Stockfish's own training notes "
+                          "report good results around 2.6 -- widens the "
+                          "loss ratio between small and large WDL-space "
+                          "errors, pushing harder to fix decisively-wrong "
+                          "predictions. See powered_loss()'s docstring.")
+    ap.add_argument("--material-loss-weight", type=float, default=0.0,
+                     help="weight on the auxiliary material_probe MSE "
+                          "term (plain pawns^2, NOT WDL-space -- so this "
+                          "is typically orders of magnitude larger in raw "
+                          "scale than the WDL loss; start small, e.g. "
+                          "0.001-0.01, and check sanity_check_net.py "
+                          "rather than assume a value works). 0.0 "
+                          "(default) disables it entirely -- the probe "
+                          "head still exists but gets no gradient, "
+                          "equivalent to not having added it.")
+    ap.add_argument("--threads", type=int, default=16,
+                     help="torch CPU thread count -- kept below this box's "
+                          "24 cores (not maxed out) since this is a shared "
+                          "multi-user box, matching the --workers default "
+                          "used elsewhere in this pipeline. Ignored when "
+                          "--device is not cpu.")
+    ap.add_argument("--device", default="cpu",
+                     help="torch device, e.g. 'cpu' (default, local "
+                          "iteration -- see this module's docstring), "
+                          "'cuda' (a rented GPU box, Phase B), or 'mps' "
+                          "(Apple Silicon). 'cuda'/'mps' with no matching "
+                          "torch/driver support available fails fast at "
+                          "startup rather than silently falling back to "
+                          "cpu.")
     ap.add_argument("--eval-every-steps", type=int, default=2000)
+    ap.add_argument("--patience", type=int, default=0,
+                     help="stop early if holdout_mse (the end-of-epoch "
+                          "value, not the more frequent mid-epoch one) "
+                          "fails to improve for this many consecutive "
+                          "epochs. 0 (default) disables early stopping -- "
+                          "always runs the full --epochs. Overfitting is "
+                          "still tracked either way (see best.json);  "
+                          "this only controls whether training keeps "
+                          "burning wall-clock time after it starts.")
     ap.add_argument("--resume", default=None,
                      help="checkpoint .pt to resume from")
     ap.add_argument("--seed", type=int, default=0)
     args = ap.parse_args()
 
-    torch.set_num_threads(args.threads)
+    if args.device.startswith("cuda") and not torch.cuda.is_available():
+        raise SystemExit(f"--device {args.device!r} requested but "
+                          f"torch.cuda.is_available() is False -- check the "
+                          f"CUDA driver/torch build on this box rather than "
+                          f"silently training on CPU")
+    if args.device == "mps" and not torch.backends.mps.is_available():
+        raise SystemExit(f"--device 'mps' requested but "
+                          f"torch.backends.mps.is_available() is False -- "
+                          f"needs Apple Silicon + a torch build with MPS "
+                          f"support, rather than silently training on CPU")
+
+    if args.device == "cpu":
+        torch.set_num_threads(args.threads)
     torch.manual_seed(args.seed)
 
     outdir = Path(args.outdir)
@@ -133,15 +227,43 @@ def main():
     holdout_data = CsrTensors(args.holdout_dir)
     print(f"  train: {len(train_data)} positions, holdout: {len(holdout_data)} positions")
 
-    model = NnueNet()
+    model = NnueNet().to(args.device)
     optimizer = torch.optim.Adam(model.parameters(), lr=args.lr)
     start_epoch = 0
     global_step = 0
+    best_holdout_mse = float("inf")
+    best_epoch = None
+    epochs_without_improvement = 0
+
+    best_path = outdir / "best.json"
+    if best_path.exists():
+        # Resuming (or re-running against an outdir that already has
+        # history): load prior best-epoch state rather than restarting
+        # the "best so far" search from scratch, which would otherwise
+        # let a later, worse epoch overwrite best.json just because
+        # this process forgot what the real best was.
+        prior_best = json.loads(best_path.read_text())
+        best_holdout_mse = prior_best["holdout_mse"]
+        best_epoch = prior_best["epoch"]
 
     if args.resume:
+        # map_location="cpu" always on load, then .to(device) below --
+        # a checkpoint saved from a CUDA run must still be resumable on
+        # a CPU-only box (and vice versa) without the load itself
+        # depending on CUDA being available.
         ckpt = torch.load(args.resume, map_location="cpu", weights_only=True)
         model.load_state_dict(ckpt["model"])
+        model.to(args.device)
         optimizer.load_state_dict(ckpt["optimizer"])
+        if args.device != "cpu":
+            # Adam's per-parameter state (exp_avg/exp_avg_sq) loads as
+            # whatever device it was saved from (map_location="cpu"
+            # above) -- move it to match the model's new device, or
+            # the first optimizer.step() throws a device-mismatch error.
+            for state in optimizer.state.values():
+                for key, val in state.items():
+                    if torch.is_tensor(val):
+                        state[key] = val.to(args.device)
         start_epoch = ckpt["epoch"] + 1
         global_step = ckpt.get("global_step", 0)
         print(f"Resumed from {args.resume} at epoch {start_epoch}")
@@ -155,40 +277,53 @@ def main():
             epoch_n = 0
             t_epoch = time.time()
 
-            for w_idx, w_off, b_idx, b_off, stm, target in iterate_batches(
+            for w_idx, w_off, b_idx, b_off, stm, target, material in iterate_batches(
                     train_data, args.batch_size, shuffle=True, rng=rng):
+                w_idx, w_off, b_idx, b_off, stm, target, material = _to_device(
+                    (w_idx, w_off, b_idx, b_off, stm, target, material), args.device)
                 optimizer.zero_grad()
-                pred = model.win_probability(w_idx, w_off, b_idx, b_off, stm, k=args.k)
-                loss = nn.functional.mse_loss(pred, target)
+                raw, material_pred = model.forward_with_material(w_idx, w_off, b_idx, b_off, stm)
+                pred = torch.sigmoid(raw / args.k)
+                wdl_loss = powered_loss(pred, target, args.loss_exponent)
+                material_loss = nn_mse_sum(material_pred, material) / len(target)
+                loss = wdl_loss + args.material_loss_weight * material_loss
                 loss.backward()
                 optimizer.step()
 
                 bs = len(target)
-                epoch_loss += loss.item() * bs
+                epoch_loss += wdl_loss.item() * bs
                 epoch_n += bs
                 global_step += 1
 
                 if global_step % args.eval_every_steps == 0:
-                    holdout_mse = evaluate(model, holdout_data, args.batch_size, args.k)
+                    holdout_wdl, holdout_material = evaluate(
+                        model, holdout_data, args.batch_size, args.k,
+                        args.loss_exponent, args.device)
                     elapsed = time.time() - t_start
                     rec = {
                         "step": global_step, "epoch": epoch,
                         "train_mse_running": epoch_loss / max(epoch_n, 1),
-                        "holdout_mse": holdout_mse,
+                        "holdout_mse": holdout_wdl,
+                        "holdout_material_mse": holdout_material,
                         "elapsed_s": round(elapsed, 1),
                         "positions_seen": epoch_n + epoch * len(train_data),
                     }
                     print(f"  step {global_step} epoch {epoch} "
                           f"train_mse={rec['train_mse_running']:.5f} "
-                          f"holdout_mse={holdout_mse:.5f} "
+                          f"holdout_mse={holdout_wdl:.5f} "
+                          f"holdout_material_mse={holdout_material:.4f} "
                           f"elapsed={elapsed:.0f}s")
                     log_f.write(json.dumps(rec) + "\n")
                     log_f.flush()
 
-            holdout_mse = evaluate(model, holdout_data, args.batch_size, args.k)
+            holdout_wdl, holdout_material = evaluate(
+                model, holdout_data, args.batch_size, args.k,
+                args.loss_exponent, args.device)
+            holdout_mse = holdout_wdl
             epoch_elapsed = time.time() - t_epoch
             print(f"Epoch {epoch} done: train_mse={epoch_loss / max(epoch_n, 1):.5f} "
-                  f"holdout_mse={holdout_mse:.5f} ({epoch_elapsed:.0f}s)")
+                  f"holdout_mse={holdout_wdl:.5f} "
+                  f"holdout_material_mse={holdout_material:.4f} ({epoch_elapsed:.0f}s)")
 
             ckpt_path = outdir / f"ckpt_epoch{epoch:03d}.pt"
             torch.save({
@@ -201,7 +336,43 @@ def main():
             }, ckpt_path)
             print(f"  wrote {ckpt_path}")
 
-    print(f"Training done in {time.time() - t_start:.0f}s")
+            # Track the best-holdout epoch regardless of --patience, so
+            # a run that does nothing special still records which
+            # epoch to actually export -- see run_nnue_pipeline.py's
+            # best_checkpoint(), which reads this instead of just
+            # taking the highest epoch number (previously: always the
+            # LAST epoch trained, overfitting or not).
+            if holdout_mse < best_holdout_mse:
+                best_holdout_mse = holdout_mse
+                best_epoch = epoch
+                epochs_without_improvement = 0
+                best_path.write_text(json.dumps({
+                    "epoch": epoch,
+                    "holdout_mse": holdout_mse,
+                    "ckpt": ckpt_path.name,
+                }))
+                print(f"  new best (holdout_mse={holdout_mse:.5f})")
+            else:
+                epochs_without_improvement += 1
+                print(f"  not an improvement over epoch {best_epoch}'s "
+                      f"holdout_mse={best_holdout_mse:.5f} "
+                      f"({epochs_without_improvement} epoch(s) without "
+                      f"improvement)")
+                if args.patience and epochs_without_improvement >= args.patience:
+                    print(f"Stopping early: no holdout improvement for "
+                          f"{epochs_without_improvement} epochs "
+                          f"(--patience {args.patience})")
+                    break
+
+    print(f"Training done in {time.time() - t_start:.0f}s "
+          f"(best epoch: {best_epoch}, holdout_mse={best_holdout_mse:.5f})")
+
+    # Completion marker for run_nnue_pipeline.py's run_stage() --
+    # NOT the final epoch's checkpoint file, which --patience early
+    # stopping may legitimately mean never gets written.
+    (outdir / "train_complete.marker").write_text(json.dumps({
+        "best_epoch": best_epoch, "best_holdout_mse": best_holdout_mse,
+    }))
 
 
 if __name__ == "__main__":
diff --git a/annotate_pgns/overfit_sanity_check.py b/annotate_pgns/overfit_sanity_check.py
new file mode 100644
index 0000000..81d0a5f
--- /dev/null
+++ b/annotate_pgns/overfit_sanity_check.py
@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+"""First-principles pipeline sanity check, 2026-10-17: can the net
+even OVERFIT a small, fixed, noiseless target -- the engine's own
+classic Eval() score on a few thousand real positions -- regressed
+directly (plain MSE on raw pawn-ish units, no WDL sigmoid, no K, no
+lambda, no quiet filter, no quantization)?
+
+Motivation: five independent head-to-head matches today (two nets,
+two time controls, with and without the lazy probe) all showed NNUE
+losing decisively to classic HCE, including at matched search depth
+where ECM tactical-suite results had suggested NNUE should do well.
+Before concluding anything about net/data quality, rule out a
+structural pipeline bug the hard way: if this architecture, this
+feature-indexing, and this training loop cannot drive error to
+~zero on a target it should be trivially able to memorize (a few
+thousand noiseless (position, scalar) pairs, same units/sign
+convention as everything else in this pipeline), something is
+broken in feature indexing, accumulator math, or the training loop
+itself -- not a "real chess judgment is hard" problem.
+
+If this DOES overfit cleanly, the pipeline is vindicated at the
+mechanical level, and today's match losses are a real net-quality
+question, not a bug hunt.
+
+Usage:
+    python3 overfit_sanity_check.py --tensor-dir /tmp/overfit_tensors \
+        --epochs 400
+"""
+import argparse
+import time
+from pathlib import Path
+
+import numpy as np
+import torch
+
+from nnue_model import NnueNet
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--tensor-dir", required=True)
+    ap.add_argument("--epochs", type=int, default=400)
+    ap.add_argument("--lr", type=float, default=1e-3)
+    ap.add_argument("--seed", type=int, default=0)
+    args = ap.parse_args()
+
+    torch.manual_seed(args.seed)
+    d = Path(args.tensor_dir)
+    # nn.EmbeddingBag's `offsets` wants one START index per bag (length
+    # N) -- the CSR files store N+1 entries (with a trailing end
+    # sentinel, per build_training_tensors.py), so the last entry
+    # must be dropped here (nnue_train.py's _gather_batch does the
+    # equivalent by only ever building N-length offset arrays).
+    white_idx = torch.from_numpy(np.load(d / "white_idx.npy").astype(np.int64))
+    white_off = torch.from_numpy(np.load(d / "white_off.npy")[:-1])
+    black_idx = torch.from_numpy(np.load(d / "black_idx.npy").astype(np.int64))
+    black_off = torch.from_numpy(np.load(d / "black_off.npy")[:-1])
+    stm = torch.from_numpy(np.load(d / "stm.npy").astype(np.int64))
+    target = torch.from_numpy(np.load(d / "target_raw.npy"))
+    n = len(target)
+    print(f"Loaded {n} (position, classic_eval) pairs -- target range "
+          f"[{target.min():.2f}, {target.max():.2f}] pawns")
+
+    model = NnueNet()
+    optimizer = torch.optim.Adam(model.parameters(), lr=args.lr)
+
+    t0 = time.time()
+    for epoch in range(args.epochs):
+        optimizer.zero_grad()
+        pred = model.forward(white_idx, white_off, black_idx, black_off, stm)
+        loss = torch.nn.functional.mse_loss(pred, target)
+        loss.backward()
+        optimizer.step()
+        if epoch % 20 == 0 or epoch == args.epochs - 1:
+            rmse_cp = (loss.item() ** 0.5) * 100
+            print(f"  epoch {epoch:4d}  train MSE (pawns^2) = {loss.item():.5f}  "
+                  f"RMSE = {rmse_cp:.1f}cp  ({time.time() - t0:.0f}s)")
+
+    with torch.no_grad():
+        pred = model.forward(white_idx, white_off, black_idx, black_off, stm)
+    final_mse = torch.nn.functional.mse_loss(pred, target).item()
+    final_rmse_cp = (final_mse ** 0.5) * 100
+    max_err_cp = (pred - target).abs().max().item() * 100
+    print(f"\nFinal: train MSE = {final_mse:.5f} pawns^2, "
+          f"RMSE = {final_rmse_cp:.1f}cp, max |error| = {max_err_cp:.1f}cp")
+    print("Verdict: ", end="")
+    if final_rmse_cp < 10:
+        print("CLEAN OVERFIT -- pipeline mechanically sound. "
+              "Today's match losses are a net-quality/data question, not a bug.")
+    elif final_rmse_cp < 50:
+        print("PARTIAL overfit -- some signal getting through, but not "
+              "cleanly memorizing a noiseless target. Worth investigating "
+              "further (LR, architecture capacity, or a real bug) before "
+              "trusting any training result built on this pipeline.")
+    else:
+        print("FAILED to overfit even a tiny noiseless target -- strong "
+              "evidence of a structural bug (feature indexing, accumulator "
+              "math, or training loop), not a net-quality/data problem. "
+              "Do not trust any prior training result until this is found.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/annotate_pgns/parse_stockfish_net.py b/annotate_pgns/parse_stockfish_net.py
new file mode 100644
index 0000000..ad70eb9
--- /dev/null
+++ b/annotate_pgns/parse_stockfish_net.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+"""Extracts the FeatureTransformer's psqtWeights table from a real
+Stockfish .nnue file (nn-fdc1d0fe6455.nnue, added to Stockfish 2023-06-12,
+rev 38e61663d836e062af0bc002814ad5149c4b7729 -- see
+https://github.com/official-stockfish/networks).
+
+Parsed against that exact commit's header-reading code (fetched and
+read directly, not guessed/reverse-engineered by eye -- see
+src/nnue/evaluate_nnue.cpp:94-125, src/nnue/nnue_feature_transformer.h:254-258,
+src/nnue/nnue_architecture.h:44, src/nnue/features/half_ka_v2_hm.h):
+
+Top-level header (evaluate_nnue.cpp read_header):
+    version   uint32 LE
+    hashValue uint32 LE
+    size      uint32 LE  (length of description string)
+    description  `size` bytes
+
+FeatureTransformer section (Detail::read_parameters + its own
+read_parameters):
+    ft_header_hash  uint32 LE, must equal
+                    HalfKAv2_hm::HashValue (0x7f234cb8) XOR (2*L1)
+    biases          L1 int16 values
+    weights         L1 * InputDimensions int16 values
+    psqtWeights     PSQTBuckets * InputDimensions int32 values,
+                    laid out psqtWeights[index*PSQTBuckets + k]
+
+InputDimensions = SQUARE_NB*PS_NB/2 = 22528, IDENTICAL formula/indexing
+to typhoon's own HalfKAv2_hm spec (docs/NNUE.md) -- same king buckets,
+same orientation table, same piece-square planes. PSQTBuckets=8
+(nnue_architecture.h). L1 for this net is NOT assumed -- derived below
+from the stored ft hash, since Stockfish's own loader treats that
+hash as the authoritative check, not file size or guessed offsets.
+
+This script ONLY extracts psqtWeights (a literal per-feature-index
+scalar linear value, one per material-count bucket) for comparison/
+diagnostic purposes against our own training coverage gaps -- NOT for
+copying into our net's weights. See the conversation: the goal is a
+from-scratch net, this is a sanity check on whether our dusty corners
+are in the right ballpark, not a transplant.
+
+Usage:
+    python3 parse_stockfish_net.py nn-fdc1d0fe6455.nnue --out sf_psqt.npy
+"""
+import argparse
+import struct
+import sys
+
+import numpy as np
+
+HALFKAV2_HM_HASH = 0x7F234CB8
+PS_NB = 11 * 64  # 704
+INPUT_DIMENSIONS = (64 * PS_NB) // 2  # 22528
+PSQT_BUCKETS = 8
+
+
+def read_u32(f):
+    return struct.unpack("<I", f.read(4))[0]
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("nnue_file")
+    ap.add_argument("--out", required=True,
+                     help="output .npy, shape [INPUT_DIMENSIONS, PSQT_BUCKETS] int32")
+    args = ap.parse_args()
+
+    with open(args.nnue_file, "rb") as f:
+        version = read_u32(f)
+        top_hash = read_u32(f)
+        desc_len = read_u32(f)
+        desc = f.read(desc_len).decode("utf-8", errors="replace")
+        print(f"version=0x{version:08x} top_hash=0x{top_hash:08x}")
+        print(f"description: {desc!r}")
+
+        ft_hash = read_u32(f)
+        print(f"ft_hash=0x{ft_hash:08x}")
+
+        # ft_hash = FeatureSet::HashValue ^ (OutputDimensions * 2)
+        # (feature_transformer.h:250) -- solve for L1 (OutputDimensions).
+        l1_times_2 = ft_hash ^ HALFKAV2_HM_HASH
+        if l1_times_2 % 2 != 0:
+            sys.exit(f"ft_hash XOR HalfKAv2_hm hash gave an odd value "
+                      f"({l1_times_2}) -- this net's feature set isn't "
+                      f"HalfKAv2_hm, or the hash constant is wrong for "
+                      f"this Stockfish version. Refusing to guess further.")
+        l1 = l1_times_2 // 2
+        print(f"Derived L1 (TransformedFeatureDimensions) = {l1}")
+
+        n_biases = l1
+        n_weights = l1 * INPUT_DIMENSIONS
+        n_psqt = PSQT_BUCKETS * INPUT_DIMENSIONS
+
+        biases = np.frombuffer(f.read(2 * n_biases), dtype="<i2")
+        assert len(biases) == n_biases, "short read on biases -- bad offset math"
+        print(f"biases: read {len(biases)}, "
+              f"range [{biases.min()}, {biases.max()}]")
+
+        # Skip the weights array -- not needed for this diagnostic, and
+        # it's ~66MB (l1 * 22528 * 2 bytes), the bulk of the file.
+        f.seek(2 * n_weights, 1)
+
+        psqt_raw = np.frombuffer(f.read(4 * n_psqt), dtype="<i4")
+        if len(psqt_raw) != n_psqt:
+            sys.exit(f"short read on psqtWeights: got {len(psqt_raw)}, "
+                      f"expected {n_psqt} -- offset math above must be "
+                      f"wrong (L1 derivation or file truncated)")
+
+    psqt = psqt_raw.reshape(INPUT_DIMENSIONS, PSQT_BUCKETS)
+    print(f"psqtWeights: shape {psqt.shape}, "
+          f"range [{psqt.min()}, {psqt.max()}]")
+    np.save(args.out, psqt)
+    print(f"Wrote {args.out}")
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/annotate_pgns/quiet_filter_games.py b/annotate_pgns/quiet_filter_games.py
index f4de31e..6b1fc58 100644
--- a/annotate_pgns/quiet_filter_games.py
+++ b/annotate_pgns/quiet_filter_games.py
@@ -24,6 +24,7 @@ Usage:
 import argparse
 import json
 import sys
+from multiprocessing import Pool
 
 from quiet_common import chess, position_is_quiet
 from quiet_filter import pv_first_move
@@ -54,10 +55,30 @@ def filter_game_positions(game_rec):
     return kept, n_excluded_check, n_excluded_capture, n_indeterminate
 
 
+def filter_line(line):
+    """Runs in a worker process: filters one game-record JSON line,
+    returning the re-serialized output line (None if the game was
+    dropped entirely) plus counts, in the same order lines came in."""
+    game_rec = json.loads(line)
+    n_pos_in = len(game_rec["positions"])
+    kept, c, p, i = filter_game_positions(game_rec)
+    if not kept:
+        return None, n_pos_in, 0, c, p, i
+    game_rec["positions"] = kept
+    return json.dumps(game_rec), n_pos_in, len(kept), c, p, i
+
+
 def main():
     ap = argparse.ArgumentParser(description=__doc__)
     ap.add_argument("--in", dest="infile", required=True)
     ap.add_argument("--out", dest="outfile", required=True)
+    ap.add_argument("--workers", type=int, default=16,
+                     help="worker processes for the per-game filtering "
+                          "pass (each game is independent); 1 disables "
+                          "multiprocessing")
+    ap.add_argument("--chunk-size", type=int, default=2000,
+                     help="games per unit of work handed to a worker; "
+                          "only relevant when --workers > 1")
     args = ap.parse_args()
 
     n_games_in = n_games_out = n_games_dropped = 0
@@ -65,24 +86,24 @@ def main():
     n_excluded_check = n_excluded_capture = n_indeterminate = 0
 
     with open(args.infile) as fin, open(args.outfile, "w") as fout:
-        for line in fin:
-            game_rec = json.loads(line)
+        if args.workers == 1:
+            results = (filter_line(line) for line in fin)
+        else:
+            with Pool(processes=args.workers) as pool:
+                results = list(pool.imap(filter_line, fin,
+                                          chunksize=args.chunk_size))
+        for out_line, pin, pout, c, p, i in results:
             n_games_in += 1
-            n_pos_in += len(game_rec["positions"])
-
-            kept, c, p, i = filter_game_positions(game_rec)
+            n_pos_in += pin
             n_excluded_check += c
             n_excluded_capture += p
             n_indeterminate += i
-
-            if not kept:
+            if out_line is None:
                 n_games_dropped += 1
                 continue
-
-            game_rec["positions"] = kept
-            fout.write(json.dumps(game_rec) + "\n")
+            fout.write(out_line + "\n")
             n_games_out += 1
-            n_pos_out += len(kept)
+            n_pos_out += pout
 
             if n_games_in % 2000 == 0:
                 print(f"  ...{n_games_in} games processed "
diff --git a/annotate_pgns/run_nnue_pipeline.py b/annotate_pgns/run_nnue_pipeline.py
new file mode 100644
index 0000000..39bc487
--- /dev/null
+++ b/annotate_pgns/run_nnue_pipeline.py
@@ -0,0 +1,382 @@
+#!/usr/bin/env python3
+"""End-to-end NNUE training pipeline driver: feed it a glob of PGN
+files that already carry per-move score comments (match_play.py's
+"{+-}NN.NN/DEPTH PV..." format -- e.g. pgn/stockfish_self_play_training_data.pgn,
+pgn/self_play_games.pgn) and (after a long time -- extract is cheap,
+but a from-scratch training run over millions of positions is not)
+it drives every remaining stage through to a .nnue weights file the
+C engine can load.
+
+This does NOT run the engine's own sd-N search (that is annotate.py's
+job, for PGNs that do NOT already have scores -- see that script and
+extract_positions.py/build_games.py instead if that is what's needed).
+Every input PGN here already has its own score/depth/PV comment on
+each move, so the only extraction step is parsing those comments --
+extract_selfplay_positions.py, despite its name, is a general reader
+for this comment format regardless of which engine(s) produced it
+(see that script's own MATE_CLAMP_THRESHOLD comment, added
+specifically to handle a foreign-engine/Stockfish-sourced corpus).
+
+Chains the existing, already-individually-tested pipeline scripts in
+this order:
+
+    1. extract_selfplay_positions.py  each input PGN -> games.jsonl
+                                       (one run per file, since that
+                                       script takes a single PGN;
+                                       outputs concatenated)
+    2. split_holdout.py               games.jsonl -> train/holdout
+                                       split (by game, not position)
+    3. fit_k_and_label.py             fits the WDL-sigmoid K on
+                                       train, labels both splits
+                                       (sharded, numpy-vectorized)
+    4. quiet_filter_games.py          drops in-check / SEE-positive-
+                                       capture positions from both
+                                       labeled splits
+    5. build_training_tensors.py      quiet-filtered games -> cached
+                                       CSR feature tensors (sharded)
+    6. nnue_train.py                  trains a checkpoint per epoch
+    7. export_nnue_weights.py         latest checkpoint -> flat
+                                       .nnue binary
+
+Each stage is skipped if its expected output already exists (pass
+--force to re-run everything regardless) -- stages 5-6 here can run
+for hours, so a crash/interrupt partway through should not mean
+starting over from stage 1.
+
+This is a thin orchestrator: it shells out to the existing, individually
+tested scripts rather than reimplementing any of their logic, and
+assumes they are run from this directory (annotate_pgns/).
+
+Usage:
+    python3 run_nnue_pipeline.py \\
+        --pgn-glob '../pgn/stockfish*.pgn' \\
+        --rundir runs/stockfish_labeled \\
+        --outfile runs/stockfish_labeled/nnue_weights.nnue
+"""
+import argparse
+import glob
+import json
+import subprocess
+import sys
+import time
+from pathlib import Path
+
+HERE = Path(__file__).resolve().parent
+
+
+def split_pgn(path, n, outdir):
+    """Splits one PGN file into n roughly-equal-game-count chunk files
+    under outdir (named <path.stem>.chunk000.pgn etc.), returning their
+    paths. Game boundaries are found by scanning for lines starting
+    "[Event " (every game in this corpus starts with one -- see
+    extract_selfplay_positions.py's own game headers,
+    build_pgn_game() in match_play.py), tracked as raw byte offsets via
+    f.buffer.tell() rather than parsed/re-rendered by python-chess --
+    avoids any risk of a PGN round-trip through its writer subtly
+    reformatting a comment or header in a way
+    extract_selfplay_positions.py's own parse wouldn't expect, and
+    avoids the text-mode tell()-vs-byte-offset mismatch a
+    chess.pgn.read_game-based scan would hit. A single-PGN-file run
+    (98k games in one file, 2026-10-19's first real pipeline attempt)
+    was otherwise stuck running extract_selfplay_positions.py
+    single-threaded for well over an hour -- that parse is purely
+    per-game, embarrassingly parallel, just never sharded because the
+    pipeline only ever sharded across *separate input files*, not
+    within one big one."""
+    offsets = []
+    with open(path, "rb") as f:
+        pos = 0
+        for raw_line in f:
+            if raw_line.startswith(b"[Event "):
+                offsets.append(pos)
+            pos += len(raw_line)
+        end = pos
+
+    n_games = len(offsets)
+    if n_games == 0:
+        return [path]
+    n = min(n, n_games)
+    boundaries = [offsets[round(i * n_games / n)] for i in range(n)] + [end]
+
+    outdir = Path(outdir)
+    outdir.mkdir(parents=True, exist_ok=True)
+    stem = Path(path).stem
+    chunk_paths = []
+    with open(path, "rb") as f:
+        for i in range(n):
+            start, stop = boundaries[i], boundaries[i + 1]
+            f.seek(start)
+            data = f.read(stop - start)
+            chunk_path = outdir / f"{stem}.chunk{i:03d}.pgn"
+            with open(chunk_path, "wb") as out:
+                out.write(data)
+            chunk_paths.append(chunk_path)
+    return chunk_paths
+
+
+def run_stage(name, outputs, cmd, force):
+    """Runs cmd (a list of str/Path args) unless every path in
+    outputs already exists and force is False. Exits the whole
+    pipeline on a non-zero return code -- no point starting a later
+    stage against a broken earlier one."""
+    outputs = [Path(p) for p in outputs]
+    if not force and all(p.exists() for p in outputs):
+        print(f"[skip] {name}: all outputs already exist "
+              f"({', '.join(str(p) for p in outputs)})")
+        return
+
+    print(f"[run]  {name}")
+    print(f"       {' '.join(str(c) for c in cmd)}")
+    t0 = time.time()
+    subprocess.run([str(c) for c in cmd], check=True)
+    print(f"[done] {name} ({time.time() - t0:.0f}s)")
+
+
+def best_checkpoint(ckpt_dir):
+    """Returns the best-holdout checkpoint (nnue_train.py's best.json),
+    not just the highest epoch number -- a run that overfits past its
+    best epoch would otherwise silently export the worst-of-the-run
+    checkpoint just because it happened to be trained last."""
+    ckpt_dir = Path(ckpt_dir)
+    best_path = ckpt_dir / "best.json"
+    if best_path.exists():
+        best = json.loads(best_path.read_text())
+        ckpt_path = ckpt_dir / best["ckpt"]
+        print(f"[info] best checkpoint: epoch {best['epoch']}, "
+              f"holdout_mse={best['holdout_mse']:.5f} ({ckpt_path})")
+        return ckpt_path
+    # Fallback for a ckpt_dir produced before best.json existed.
+    ckpts = sorted(ckpt_dir.glob("ckpt_epoch*.pt"))
+    if not ckpts:
+        sys.exit(f"no checkpoints found in {ckpt_dir} -- nnue_train.py "
+                  f"stage did not produce any?")
+    print(f"[warn] no best.json in {ckpt_dir} -- falling back to the "
+          f"last-trained epoch, not necessarily the best holdout one")
+    return ckpts[-1]
+
+
+def main():
+    ap = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--pgn-glob", action="append", required=True,
+                     help="glob pattern (quoted, so this script expands "
+                          "it rather than the shell) for input PGN "
+                          "files, each already carrying per-move "
+                          "score/depth/PV comments; repeatable")
+    ap.add_argument("--rundir", required=True,
+                     help="directory for all intermediate/output "
+                          "artifacts of this run")
+    ap.add_argument("--outfile", required=True,
+                     help="final .nnue weights file path")
+    ap.add_argument("--event", default=None,
+                     help="extract_selfplay_positions.py: only extract "
+                          "games with this Event header (default: any)")
+    ap.add_argument("--dates", nargs="*", default=None,
+                     help="extract_selfplay_positions.py: only extract "
+                          "games with Date header in this list "
+                          "(default: any)")
+    ap.add_argument("--exclude-movers", nargs="*", default=[],
+                     help="extract_selfplay_positions.py: skip plies "
+                          "whose mover's PGN name is in this list")
+    ap.add_argument("--holdout-frac", type=float, default=0.1,
+                     help="split_holdout.py's game-level holdout "
+                          "fraction")
+    ap.add_argument("--workers", type=int, default=16,
+                     help="worker processes for every sharded Python "
+                          "stage (fit_k_and_label.py, "
+                          "build_training_tensors.py). Default 16, "
+                          "kept below this box's 24 cores to leave "
+                          "other jobs some breathing room")
+    ap.add_argument("--threads", type=int, default=16,
+                     help="torch CPU thread count for nnue_train.py, "
+                          "same reasoning as --workers (ignored if "
+                          "--device is not cpu)")
+    ap.add_argument("--device", default="cpu",
+                     help="nnue_train.py's --device, e.g. 'cuda' on a "
+                          "rented GPU box")
+    ap.add_argument("--lambda", dest="lam", type=float, default=0.9,
+                     help="build_training_tensors.py's "
+                          "target = lambda*p_search + (1-lambda)*r_outcome "
+                          "for non-exact rows")
+    ap.add_argument("--epochs", type=int, default=20)
+    ap.add_argument("--patience", type=int, default=0,
+                     help="nnue_train.py's --patience (0 disables early "
+                          "stopping, always runs the full --epochs)")
+    ap.add_argument("--batch-size", type=int, default=8192)
+    ap.add_argument("--lr", type=float, default=1e-3)
+    ap.add_argument("--loss-exponent", type=float, default=2.0)
+    ap.add_argument("--material-loss-weight", type=float, default=0.0)
+    ap.add_argument("--force", action="store_true",
+                     help="re-run every stage even if its outputs "
+                          "already exist")
+    args = ap.parse_args()
+
+    pgn_paths = sorted({p for pattern in args.pgn_glob
+                         for p in glob.glob(pattern)})
+    if not pgn_paths:
+        sys.exit(f"--pgn-glob matched no files: {args.pgn_glob!r}")
+    print(f"Input PGNs ({len(pgn_paths)}): {pgn_paths}")
+
+    rundir = Path(args.rundir)
+    rundir.mkdir(parents=True, exist_ok=True)
+
+    games_path = rundir / "games.jsonl"
+    train_path = rundir / "games.train.jsonl"
+    holdout_path = rundir / "games.holdout.jsonl"
+    train_labeled_path = rundir / "games.train.labeled.jsonl"
+    holdout_labeled_path = rundir / "games.holdout.labeled.jsonl"
+    train_quiet_path = rundir / "games.train.labeled.quiet.jsonl"
+    holdout_quiet_path = rundir / "games.holdout.labeled.quiet.jsonl"
+    tensors_train_dir = rundir / "tensors_train"
+    tensors_holdout_dir = rundir / "tensors_holdout"
+    ckpt_dir = rundir / "nnue_ckpt"
+
+    event_flag = ["--event", args.event] if args.event is not None else []
+    dates_flag = ["--dates", *args.dates] if args.dates else []
+    exclude_flag = (["--exclude-movers", *args.exclude_movers]
+                     if args.exclude_movers else [])
+
+    # 1. Extract already-scored positions. extract_selfplay_positions.py
+    # is purely per-game, embarrassingly parallel -- shard each input
+    # PGN into args.workers game-aligned chunks (split_pgn, raw byte
+    # copy, no PGN re-render) and run one extraction process per chunk
+    # concurrently, rather than one single-threaded process per whole
+    # input file (which is what left a 98k-game single-file run stuck
+    # for over an hour on 2026-10-19).
+    split_dir = rundir / "pgn_splits"
+    all_chunks = []  # (chunk_path, part_path) pairs, in output order
+    for pgn_path in pgn_paths:
+        chunks = split_pgn(pgn_path, args.workers, split_dir)
+        for chunk_path in chunks:
+            part_path = rundir / f"{chunk_path.stem}.jsonl"
+            all_chunks.append((chunk_path, part_path))
+
+    pending = [(c, p) for c, p in all_chunks if args.force or not p.exists()]
+    if pending:
+        print(f"[run]  extract_selfplay_positions.py: {len(pending)} chunk(s) "
+              f"across {len(pgn_paths)} input file(s), up to {args.workers} "
+              f"concurrent")
+        t0 = time.time()
+        procs = []
+        for chunk_path, part_path in pending:
+            cmd = [sys.executable, HERE / "extract_selfplay_positions.py",
+                   chunk_path, "--out", part_path,
+                   *event_flag, *dates_flag, *exclude_flag]
+            while len(procs) >= args.workers:
+                procs = [p for p in procs if p.poll() is None]
+                if len(procs) >= args.workers:
+                    time.sleep(0.5)
+            procs.append(subprocess.Popen([str(c) for c in cmd]))
+        for p in procs:
+            p.wait()
+        failed = [p for p in procs if p.returncode != 0]
+        if failed:
+            sys.exit(f"extract_selfplay_positions.py: {len(failed)} chunk(s) "
+                      f"failed -- aborting pipeline")
+        print(f"[done] extract_selfplay_positions.py ({time.time() - t0:.0f}s)")
+    else:
+        print(f"[skip] extract_selfplay_positions.py: all "
+              f"{len(all_chunks)} chunk outputs already exist")
+
+    if args.force or not games_path.exists():
+        print(f"[run]  concatenating {len(all_chunks)} parts -> {games_path}")
+        with open(games_path, "w") as out:
+            for _, part_path in all_chunks:
+                with open(part_path) as f:
+                    out.write(f.read())
+    else:
+        print(f"[skip] concatenate: {games_path} already exists")
+
+    # 2. Holdout split (by game).
+    run_stage(
+        "split_holdout.py", [train_path, holdout_path],
+        [sys.executable, HERE / "split_holdout.py",
+         "--in", games_path, "--train-out", train_path,
+         "--holdout-out", holdout_path,
+         "--holdout-frac", args.holdout_frac],
+        args.force)
+
+    # 3. Fit K + label both splits (sharded, numpy-vectorized).
+    run_stage(
+        "fit_k_and_label.py", [train_labeled_path, holdout_labeled_path],
+        [sys.executable, HERE / "fit_k_and_label.py",
+         "--train", train_path, "--holdout", holdout_path,
+         "--train-out", train_labeled_path,
+         "--holdout-out", holdout_labeled_path,
+         "--workers", args.workers],
+        args.force)
+
+    # 4. Quiet-filter both labeled splits.
+    run_stage(
+        "quiet_filter_games.py (train)", [train_quiet_path],
+        [sys.executable, HERE / "quiet_filter_games.py",
+         "--in", train_labeled_path, "--out", train_quiet_path,
+         "--workers", args.workers],
+        args.force)
+    run_stage(
+        "quiet_filter_games.py (holdout)", [holdout_quiet_path],
+        [sys.executable, HERE / "quiet_filter_games.py",
+         "--in", holdout_labeled_path, "--out", holdout_quiet_path,
+         "--workers", args.workers],
+        args.force)
+
+    # 5. Build cached CSR feature tensors (sharded).
+    run_stage(
+        "build_training_tensors.py (train)",
+        [tensors_train_dir / "target.npy"],
+        [sys.executable, HERE / "build_training_tensors.py",
+         "--in", train_quiet_path, "--outdir", tensors_train_dir,
+         "--lambda", args.lam, "--workers", args.workers],
+        args.force)
+    run_stage(
+        "build_training_tensors.py (holdout)",
+        [tensors_holdout_dir / "target.npy"],
+        [sys.executable, HERE / "build_training_tensors.py",
+         "--in", holdout_quiet_path, "--outdir", tensors_holdout_dir,
+         "--lambda", args.lam, "--workers", args.workers],
+        args.force)
+
+    # 6. Train. Marker: train_complete.marker, written by nnue_train.py
+    # itself right at the end of main() -- NOT the final epoch's
+    # checkpoint file, which --patience early stopping may legitimately
+    # mean never gets written, and which would otherwise make a clean
+    # early-stopped run look unfinished forever to this existence
+    # check. If a previous attempt left partial epoch checkpoints
+    # without that marker (killed partway, or crashed), resume from the
+    # latest one instead of silently restarting at epoch 0 and
+    # discarding completed epochs.
+    train_complete_marker = ckpt_dir / "train_complete.marker"
+    resume_flag = []
+    if not args.force and not train_complete_marker.exists():
+        existing_ckpts = sorted(ckpt_dir.glob("ckpt_epoch*.pt"))
+        if existing_ckpts:
+            resume_flag = ["--resume", existing_ckpts[-1]]
+            print(f"[info] nnue_train.py: resuming from "
+                  f"{existing_ckpts[-1]} (no train_complete.marker found)")
+    run_stage(
+        "nnue_train.py", [train_complete_marker],
+        [sys.executable, HERE / "nnue_train.py",
+         "--train-dir", tensors_train_dir,
+         "--holdout-dir", tensors_holdout_dir,
+         "--outdir", ckpt_dir, *resume_flag,
+         "--epochs", args.epochs, "--batch-size", args.batch_size,
+         "--lr", args.lr, "--loss-exponent", args.loss_exponent,
+         "--material-loss-weight", args.material_loss_weight,
+         "--threads", args.threads, "--device", args.device,
+         "--patience", args.patience],
+        args.force)
+
+    # 7. Export the latest checkpoint to the flat .nnue binary format.
+    ckpt = best_checkpoint(ckpt_dir)
+    run_stage(
+        "export_nnue_weights.py", [args.outfile],
+        [sys.executable, HERE / "export_nnue_weights.py",
+         ckpt, args.outfile, "--tensor-dir", tensors_train_dir],
+        args.force)
+
+    print(f"\nPipeline complete: {args.outfile}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/annotate_pgns/sample_egtb_positions.py b/annotate_pgns/sample_egtb_positions.py
new file mode 100644
index 0000000..f4cfa45
--- /dev/null
+++ b/annotate_pgns/sample_egtb_positions.py
@@ -0,0 +1,202 @@
+#!/usr/bin/env python3
+"""Generate direct Syzygy-ground-truth training positions for the
+EGTB-mixing idea in NNUE.md: random LEGAL positions for a curated set
+of endgame piece configurations, each probed directly against the
+Syzygy tablebase for its exact WDL result and labeled with that result
+(1.0/0.0/0.5, mover-relative) -- NOT derived from any engine search
+score, so there's no sigmoid/lambda involved at all for these rows
+(see build_training_tensors.py's position_target: exact rows with an
+r_outcome use it directly).
+
+Why this exists as a separate generator rather than just trusting
+Stockfish self-play games to supply this signal: a self-play game's
+score field, even when the search is quietly influenced by a TB probe
+a few plies down, is NOT a reliable way to detect "this exact position
+is a TB-determined result" after the fact -- see parse_comment's
+MATE_CLAMP_THRESHOLD handling in extract_selfplay_positions.py, which
+can only catch the unambiguous extreme case (|score| close to the
+engine's own mate-score cap). An ordinary-looking, e.g., "+0.30"
+6-piece position could just as easily be a genuine TB draw as a normal
+balanced middlegame-ish score, and there is no way to tell the two
+apart from the log alone. Rather than guess, this script sidesteps the
+ambiguity entirely: it builds positions whose piece count and
+side-to-move are known up front, and asks the tablebase directly,
+which is unambiguous by construction.
+
+Curated piece-set list: endgames that are either common in real games
+(KPK, KRK, KQK, basic mating/drawing technique) or specifically
+instructive for a net's positional judgment near the piece-trade
+boundary (KRPKR, KBPKB, KNPKN -- same-minor-vs-pawn fortress/win
+distinctions; KPPKPP -- connected/outside-passer races). Deliberately
+excludes same-type-heavy, essentially-never-arising-from-real-play
+configurations like KNNNKB (three knights) or KQQKQQ -- these exist in
+the 6-piece Syzygy set for completeness but add no realistic training
+signal, only wasted sampling effort chasing positions the net will
+never see at inference time.
+
+Usage:
+    python3 sample_egtb_positions.py \\
+        --egtb-path /zscratch/egtb \\
+        --out runs/egtb_direct/positions.jsonl \\
+        --per-config 2000
+"""
+import argparse
+import json
+import random
+import sys
+
+import chess
+import chess.syzygy
+
+# Each entry: (white_pieces, black_pieces), piece letters excluding the
+# king (every side always has exactly one king, added automatically).
+# Ordered white-majority-first by convention; sampling itself places
+# White/Black and side-to-move uniformly at random regardless.
+USEFUL_CONFIGS = [
+    # 3-piece: the universal basics.
+    ("P", ""),      # KPK
+    ("R", ""),      # KRK
+    ("Q", ""),      # KQK
+    ("B", "B"),     # KBKB (same-color-bishop draw / opposite-color technique)
+    ("N", "N"),     # KNKN
+    # 4-piece: one side up material, or pawn-down defensive technique.
+    ("R", "P"),     # KRKP
+    ("Q", "P"),     # KQKP
+    ("Q", "R"),     # KQKR
+    ("R", "B"),     # KRKB
+    ("R", "N"),     # KRKN
+    ("B", "N"),     # KBKN
+    ("P", "P"),     # KPKP (opposition/zugzwang pawn races)
+    # 5-piece: minor+pawn fortress/technique boundary cases, specifically
+    # called out as instructive in the discussion that motivated this.
+    ("RP", "R"),    # KRPKR (the single most common rook-endgame skeleton)
+    ("BP", "B"),    # KBPKB (same/opposite-color bishop-vs-pawn technique)
+    ("NP", "N"),    # KNPKN
+    ("QP", "Q"),    # KQPKQ
+    ("PP", "P"),    # KPPKP
+    # 6-piece: the specific "KPPKPP, not KNNNKB" case from the discussion.
+    ("PP", "PP"),   # KPPKPP (connected/outside passer races)
+    ("RP", "RP"),   # KRPKRP
+    ("BP", "BP"),   # KBPKBP
+    ("RR", "R"),    # KRRKR
+]
+
+PIECE_SYMBOL = {
+    "P": chess.PAWN, "N": chess.KNIGHT, "B": chess.BISHOP,
+    "R": chess.ROOK, "Q": chess.QUEEN,
+}
+
+
+def random_board(white_pieces, black_pieces, rng):
+    """Attempts one random legal placement for this piece configuration
+    (kings included). Returns a chess.Board, or None if this attempt
+    didn't produce a legal position (retried by the caller) -- most
+    retries are from kings ending up adjacent or the side not on move
+    being left in check, both of which are just rejection-sampled away
+    rather than avoided by construction, since legal-position rejection
+    sampling is simple, fast for this few-piece count, and impossible
+    to get subtly wrong the way a hand-built adjacency/legality check
+    could."""
+    squares = rng.sample(range(64), 2 + len(white_pieces) + len(black_pieces))
+    board = chess.Board.empty()
+    board.set_piece_at(squares[0], chess.Piece(chess.KING, chess.WHITE))
+    board.set_piece_at(squares[1], chess.Piece(chess.KING, chess.BLACK))
+    idx = 2
+    for letter in white_pieces:
+        sq = squares[idx]
+        idx += 1
+        if letter == "P" and chess.square_rank(sq) in (0, 7):
+            return None
+        board.set_piece_at(sq, chess.Piece(PIECE_SYMBOL[letter], chess.WHITE))
+    for letter in black_pieces:
+        sq = squares[idx]
+        idx += 1
+        if letter == "P" and chess.square_rank(sq) in (0, 7):
+            return None
+        board.set_piece_at(sq, chess.Piece(PIECE_SYMBOL[letter], chess.BLACK))
+    board.turn = chess.WHITE if rng.random() < 0.5 else chess.BLACK
+    board.castling_rights = chess.BB_EMPTY
+    if not board.is_valid():
+        return None
+    return board
+
+
+def sample_config(tablebase, white_pieces, black_pieces, n, rng, max_attempts):
+    """Returns up to n (fen, r_outcome) pairs for one piece configuration,
+    r_outcome mover-relative (1.0/0.5/0.0), cursed-win/blessed-loss
+    collapsed to the practical 1.0/0.0 (not 0.5) -- under the engine's
+    own 50-move-rule handling a cursed win is still winnable with
+    correct play short of the counter running out, so it is a much
+    better training label as "winning" than as "drawn"; collapsing it
+    the other way would teach the net to misjudge a real winning
+    chance as hopeless."""
+    out = []
+    attempts = 0
+    while len(out) < n and attempts < max_attempts:
+        attempts += 1
+        board = random_board(white_pieces, black_pieces, rng)
+        if board is None:
+            continue
+        try:
+            wdl = tablebase.probe_wdl(board)
+        except chess.syzygy.MissingTableError:
+            return out  # this exact configuration has no table -- give up
+        if wdl > 0:
+            r_outcome = 1.0
+        elif wdl < 0:
+            r_outcome = 0.0
+        else:
+            r_outcome = 0.5
+        out.append((board.fen(), r_outcome))
+    return out
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--egtb-path", default="/zscratch/egtb")
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--per-config", type=int, default=2000,
+                     help="target sample count per piece configuration "
+                          "(each configuration is also sampled with both "
+                          "colors' roles swapped, so total yield per entry "
+                          "in USEFUL_CONFIGS is up to 2x this)")
+    ap.add_argument("--max-attempts-factor", type=int, default=50,
+                     help="give up on a configuration after "
+                          "per-config * this many random placements, in "
+                          "case legal-position rejection sampling is "
+                          "pathologically slow for some piece set")
+    ap.add_argument("--seed", type=int, default=0)
+    args = ap.parse_args()
+
+    rng = random.Random(args.seed)
+    tablebase = chess.syzygy.open_tablebase(args.egtb_path)
+
+    total = 0
+    with open(args.out, "w") as fout:
+        for white_pieces, black_pieces in USEFUL_CONFIGS:
+            for wp, bp in ({(white_pieces, black_pieces),
+                             (black_pieces, white_pieces)}):
+                rows = sample_config(tablebase, wp, bp, args.per_config, rng,
+                                      args.per_config * args.max_attempts_factor)
+                for fen, r_outcome in rows:
+                    game_rec = {
+                        "game_id": f"egtb_direct_K{wp}K{bp}_{total}",
+                        "result": "*",
+                        "positions": [{
+                            "fen": fen,
+                            "exact": True,
+                            "exact_reason": "egtb_direct",
+                            "r_outcome": r_outcome,
+                        }],
+                    }
+                    fout.write(json.dumps(game_rec) + "\n")
+                    total += 1
+                print(f"  K{wp}vK{bp}: {len(rows)}/{args.per_config} sampled",
+                      file=sys.stderr)
+
+    print(f"Wrote {args.out}: {total} direct-EGTB positions", file=sys.stderr)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/annotate_pgns/sanity_check_net.py b/annotate_pgns/sanity_check_net.py
new file mode 100755
index 0000000..fd3b240
--- /dev/null
+++ b/annotate_pgns/sanity_check_net.py
@@ -0,0 +1,306 @@
+#!/usr/bin/env python3
+"""Material-sensitivity sanity check for a trained NNUE checkpoint.
+
+Loads a raw PyTorch checkpoint (nnue_train.py's ckpt_epochNNN.pt) and
+runs it, via nnue_model.py's forward() directly -- no export/C
+round-trip needed, so this is the earliest point in the pipeline a
+checkpoint can be checked -- against a battery of hand-picked base
+positions, each with one piece removed at a time. Reports the raw
+score delta from the base position for each removal and flags
+whether the deltas are monotonic in material value (pawn < knight ~
+bishop < rook < queen), which a net with a real concept of material
+should satisfy even though nothing forces it to structurally.
+
+This exists because of a real, measured failure found 2026-10-17:
+removing a queen, a rook, a knight, and a pawn from the same base
+position produced nearly indistinguishable score drops (-0.56 to
+-0.70) -- directionally correct (material loss always hurts) but
+with no real sense of HOW MUCH. Likely root cause: training loss is
+MSE against sigmoid(raw/K), and sigmoid's gradient vanishes once a
+position is already clearly lost/won, so a decisively-bad position
+gets little pressure to become MORE decisively bad regardless of
+how much additional material disappears. This script exists to make
+that failure mode visible at a glance on any future checkpoint,
+rather than rediscovering it by hand each time.
+
+Base positions deliberately avoid the opening setup: deleting a
+queen from the starting position is itself wildly out-of-distribution
+(no real game ever reaches that exact shape), so a net's response to
+it isn't a trustworthy read -- these use realistic midgame/endgame-ish
+material counts instead.
+
+2026-10-17: added --cross-sectional, a population-level version of the
+same idea. The hand-picked BASE_POSITIONS checks above are useful as a
+fast smoke test, but a single removal tests one specific piece on one
+specific square in one specific position -- its result is confounded
+by that piece's concrete role there (an active bishop vs. a boxed-in
+rook can legitimately score differently for positional reasons having
+nothing to do with the net's general sense of material), as observed
+directly: bishop-removal outscored rook-removal in one hand-picked
+position and the reverse held in another. --cross-sectional instead
+samples many REAL positions from a holdout pool, removes one piece of
+each type per position (when present), and pools deltas by piece
+type across the whole sample -- so per-position context should average
+out, leaving a cleaner read on the net's actual learned material
+prior (mean/median/stdev per piece type, not a single number).
+
+Usage:
+    python3 sanity_check_net.py --ckpt runs/pooled_gen2/nnue_ckpt/ckpt_epoch032.pt
+    python3 sanity_check_net.py --ckpt ... --cross-sectional \
+        --sample-file runs/pooled_gen2/games.holdout.labeled.quiet.jsonl \
+        --n-positions 2000
+"""
+import argparse
+import glob
+import json
+import random
+
+import chess
+import numpy as np
+import torch
+
+from nnue_feature_spec import active_feature_indices
+from nnue_model import NnueNet
+
+# (piece_type, approximate pawn value) -- value only used to order
+# the monotonicity check, not fed to the net.
+PIECE_VALUES = [
+    (chess.PAWN, 1),
+    (chess.KNIGHT, 3),
+    (chess.BISHOP, 3),
+    (chess.ROOK, 5),
+    (chess.QUEEN, 9),
+]
+
+BASE_POSITIONS = {
+    "midgame_balanced":
+        "r2qk2r/ppp2ppp/2n1bn2/2b1p3/4P3/2N1BN2/PPP2PPP/R2QK2R w KQkq - 0 1",
+    "midgame_balanced_black_to_move":
+        "r2qk2r/ppp2ppp/2n1bn2/2b1p3/4P3/2N1BN2/PPP2PPP/R2QK2R b KQkq - 0 1",
+    "endgame_rook_pawns":
+        "8/5pk1/6p1/8/4R3/5PK1/6P1/8 w - - 0 1",
+    "queen_endgame":
+        "6k1/5ppp/8/8/8/8/4QPPP/6K1 w - - 0 1",
+}
+
+
+def load_model(ckpt_path):
+    ck = torch.load(ckpt_path, map_location="cpu", weights_only=True)
+    model = NnueNet()
+    model.load_state_dict(ck["model"])
+    model.eval()
+    return model, ck.get("holdout_mse")
+
+
+def score(model, board):
+    w = active_feature_indices(board, chess.WHITE)
+    b = active_feature_indices(board, chess.BLACK)
+    w_idx = torch.tensor(w, dtype=torch.long)
+    w_off = torch.tensor([0], dtype=torch.long)
+    b_idx = torch.tensor(b, dtype=torch.long)
+    b_off = torch.tensor([0], dtype=torch.long)
+    stm = torch.tensor([1 if board.turn == chess.WHITE else 0], dtype=torch.long)
+    with torch.no_grad():
+        return model.forward(w_idx, w_off, b_idx, b_off, stm).item()
+
+
+def one_sided_removals(board, color):
+    """Yields (piece_type, value, square, board_with_piece_removed)
+    for every piece of `color` on the board, one at a time -- skips
+    kings (removing a king produces an illegal/meaningless position)."""
+    for piece_type, value in PIECE_VALUES:
+        squares = list(board.pieces(piece_type, color))
+        if not squares:
+            continue
+        sq = squares[0]
+        b2 = board.copy()
+        b2.remove_piece_at(sq)
+        yield piece_type, value, sq, b2
+
+
+def run_base_position(model, name, fen):
+    board = chess.Board(fen)
+    mover_color = board.turn
+    mover_name = "White" if mover_color == chess.WHITE else "Black"
+    base_score = score(model, board)
+    print(f"\n=== {name} ({mover_name} to move) ===")
+    print(f"  base score: {base_score:+.3f}")
+
+    print(f"  -- removing {mover_name}'s OWN pieces (score should drop, "
+          f"and drop MORE for more valuable pieces) --")
+    rows = []
+    for piece_type, value, sq, b2 in one_sided_removals(board, mover_color):
+        s = score(model, b2)
+        delta = s - base_score
+        rows.append((value, chess.piece_name(piece_type), delta))
+        print(f"    remove own {chess.piece_name(piece_type):8s} "
+              f"(val={value}): score={s:+.3f}  delta={delta:+.3f}")
+
+    # Monotonicity check: sort by piece value, deltas should be
+    # non-increasing (more valuable piece removed -> more negative
+    # delta), allowing ties between knight/bishop (same nominal value).
+    rows.sort(key=lambda r: r[0])
+    deltas = [d for _, _, d in rows]
+    monotonic = all(deltas[i] >= deltas[i + 1] - 1e-6 for i in range(len(deltas) - 1))
+    spread = max(deltas) - min(deltas) if deltas else 0.0
+    print(f"  monotonic (pawn..queen delta non-increasing)? {monotonic}  "
+          f"(spread between least/most valuable piece's delta: {spread:.3f})")
+    if not monotonic:
+        print("  *** FAIL: removing a MORE valuable piece did not hurt "
+              "the score MORE -- net may not have learned material. ***")
+    if deltas and spread < 0.5:
+        print(f"  *** WARNING: spread of only {spread:.3f} across pawn..queen "
+              f"removal is small -- net may be material-insensitive even if "
+              f"technically monotonic. ***")
+
+    print(f"  -- removing opponent's pieces (score should RISE) --")
+    opp_color = not mover_color
+    for piece_type, value, sq, b2 in one_sided_removals(board, opp_color):
+        s = score(model, b2)
+        delta = s - base_score
+        flag = "" if delta > 0 else "  *** FAIL: should have risen ***"
+        print(f"    remove opp {chess.piece_name(piece_type):8s} "
+              f"(val={value}): score={s:+.3f}  delta={delta:+.3f}{flag}")
+
+
+def sample_fens(sample_file, n_positions, seed):
+    """Reservoir-sample n_positions FENs out of every position in
+    every game in a games.jsonl-shaped file (quiet_filter_games.py
+    output). Loads all FENs into memory first (a string per position,
+    cheap even at a few hundred thousand rows) rather than true
+    streaming reservoir sampling -- simpler, and the holdout pools
+    this is meant to run against are a few hundred MB at most."""
+    fens = []
+    with open(sample_file) as f:
+        for line in f:
+            game = json.loads(line)
+            for pos in game["positions"]:
+                fens.append(pos["fen"])
+    rng = random.Random(seed)
+    if n_positions < len(fens):
+        fens = rng.sample(fens, n_positions)
+    return fens
+
+
+def one_random_removal_per_piece_type(board, color, rng):
+    """Yields (piece_type, value, board_with_one_piece_removed) --
+    exactly one randomly-chosen piece per type present for `color`
+    (not every piece, unlike one_sided_removals -- across thousands
+    of sampled positions this still covers many squares per type in
+    aggregate without each position contributing >1 sample per type,
+    which would overweight positions with e.g. two bishops)."""
+    for piece_type, value in PIECE_VALUES:
+        squares = list(board.pieces(piece_type, color))
+        if not squares:
+            continue
+        sq = rng.choice(squares)
+        b2 = board.copy()
+        b2.remove_piece_at(sq)
+        yield piece_type, value, b2
+
+
+def run_cross_sectional(model, sample_file, n_positions, seed):
+    print(f"\n=== cross-sectional population check ===")
+    print(f"  sampling {n_positions} positions from {sample_file} ...")
+    fens = sample_fens(sample_file, n_positions, seed)
+    print(f"  got {len(fens)} FENs")
+
+    rng = random.Random(seed)
+    # own_deltas[piece_type] = list of deltas from removing the MOVER's
+    # own piece of that type (expected negative, more negative for
+    # more valuable pieces). opp_deltas = removing the opponent's
+    # piece instead (expected positive).
+    own_deltas = {pt: [] for pt, _ in PIECE_VALUES}
+    opp_deltas = {pt: [] for pt, _ in PIECE_VALUES}
+    n_skipped = 0
+
+    for i, fen in enumerate(fens):
+        try:
+            board = chess.Board(fen)
+        except ValueError:
+            n_skipped += 1
+            continue
+        mover = board.turn
+        base_score = score(model, board)
+        for piece_type, value, b2 in one_random_removal_per_piece_type(board, mover, rng):
+            own_deltas[piece_type].append(score(model, b2) - base_score)
+        for piece_type, value, b2 in one_random_removal_per_piece_type(board, not mover, rng):
+            opp_deltas[piece_type].append(score(model, b2) - base_score)
+        if (i + 1) % 500 == 0:
+            print(f"  ...{i + 1}/{len(fens)} positions processed")
+
+    if n_skipped:
+        print(f"  ({n_skipped} FENs failed to parse, skipped)")
+
+    def report(deltas_by_type, label, expect_sign):
+        print(f"\n  -- {label} (expect sign: {'negative' if expect_sign < 0 else 'positive'}) --")
+        means = []
+        for piece_type, value in PIECE_VALUES:
+            d = np.array(deltas_by_type[piece_type])
+            if len(d) == 0:
+                print(f"    {chess.piece_name(piece_type):8s} (val={value}): no samples")
+                continue
+            means.append((value, chess.piece_name(piece_type), d.mean()))
+            print(f"    {chess.piece_name(piece_type):8s} (val={value}): "
+                  f"n={len(d):6d}  mean={d.mean():+.4f}  median={np.median(d):+.4f}  "
+                  f"std={d.std():.4f}  p25={np.percentile(d, 25):+.4f}  "
+                  f"p75={np.percentile(d, 75):+.4f}")
+        means.sort(key=lambda r: r[0])
+        mean_vals = [m for _, _, m in means]
+        if expect_sign < 0:
+            monotonic = all(mean_vals[i] >= mean_vals[i + 1] - 1e-6 for i in range(len(mean_vals) - 1))
+        else:
+            monotonic = all(mean_vals[i] <= mean_vals[i + 1] + 1e-6 for i in range(len(mean_vals) - 1))
+        print(f"    population means monotonic in piece value? {monotonic}")
+        if not monotonic:
+            order = ", ".join(f"{name}={m:+.4f}" for _, name, m in means)
+            print(f"    *** FAIL: population-level mean deltas not ordered "
+                  f"by piece value: {order} ***")
+
+    report(own_deltas, "removing MOVER's own piece", expect_sign=-1)
+    report(opp_deltas, "removing OPPONENT's piece", expect_sign=+1)
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                  formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--ckpt", default=None,
+                     help="checkpoint path. Default: latest "
+                          "ckpt_epoch*.pt under --ckpt-glob-dir")
+    ap.add_argument("--ckpt-glob-dir", default="runs/pooled_gen2/nnue_ckpt",
+                     help="directory to find the latest checkpoint in, "
+                          "when --ckpt is not given")
+    ap.add_argument("--cross-sectional", action="store_true",
+                     help="run the population-level check (see module "
+                          "docstring's 2026-10-17 note) instead of the "
+                          "hand-picked single-position checks")
+    ap.add_argument("--sample-file",
+                     default="runs/pooled_gen2/games.holdout.labeled.quiet.jsonl",
+                     help="games.jsonl-shaped file to sample real "
+                          "positions from, for --cross-sectional")
+    ap.add_argument("--n-positions", type=int, default=2000,
+                     help="number of positions to sample, for --cross-sectional")
+    ap.add_argument("--seed", type=int, default=0)
+    args = ap.parse_args()
+
+    ckpt_path = args.ckpt
+    if ckpt_path is None:
+        candidates = sorted(glob.glob(f"{args.ckpt_glob_dir}/ckpt_epoch*.pt"))
+        if not candidates:
+            raise SystemExit(f"No checkpoints found under {args.ckpt_glob_dir}")
+        ckpt_path = candidates[-1]
+
+    print(f"Loading {ckpt_path} ...")
+    model, holdout_mse = load_model(ckpt_path)
+    if holdout_mse is not None:
+        print(f"  (checkpoint's recorded holdout_mse: {holdout_mse:.5f})")
+
+    if args.cross_sectional:
+        run_cross_sectional(model, args.sample_file, args.n_positions, args.seed)
+    else:
+        for name, fen in BASE_POSITIONS.items():
+            run_base_position(model, name, fen)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/annotate_pgns/shard_quiet_filter.py b/annotate_pgns/shard_quiet_filter.py
new file mode 100644
index 0000000..fb1c8cc
--- /dev/null
+++ b/annotate_pgns/shard_quiet_filter.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+"""Shard quiet_filter.py's input across N worker subprocesses and
+recombine their outputs -- quiet_filter.py itself is single-threaded
+(a plain line-by-line loop doing chess.Board() + SEE per row), and
+with 2.5M+ rows that's the one CPU-bound step in the post-annotation
+pipeline with no parallelism of its own, unlike build_training_tensors.py's
+--workers Pool right next to it.
+
+Deliberately wraps quiet_filter.py as a subprocess rather than
+reimplementing its filter logic here -- one definition of the filter,
+reused unmodified; this script only owns the sharding/recombining.
+
+Splits the input round-robin (line i -> shard i % --workers) in one
+pass, so no second pass to count lines first is needed. Row order in
+the combined output does NOT match the input's order (shard outputs
+are concatenated shard-by-shard) -- fine for this file's downstream
+use (positions.quiet.jsonl is a flat, order-independent training-set
+input; quiet_filter_games.py, which does care about preserving a
+game's own position order, filters the game-shaped file directly and
+does not go through this wrapper).
+
+Usage:
+    python3 shard_quiet_filter.py --in runs/sd10_twic_2M/positions.jsonl \
+        --out runs/sd10_twic_2M/positions.quiet.jsonl --workers 16
+"""
+import argparse
+import subprocess
+import sys
+import tempfile
+from pathlib import Path
+
+_QUIET_FILTER = Path(__file__).resolve().parent / "quiet_filter.py"
+
+
+def shard_input(infile, shard_paths):
+    handles = [open(p, "w") for p in shard_paths]
+    n = len(handles)
+    try:
+        with open(infile) as fin:
+            for i, line in enumerate(fin):
+                handles[i % n].write(line)
+    finally:
+        for h in handles:
+            h.close()
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--in", dest="infile", required=True)
+    ap.add_argument("--out", dest="outfile", required=True)
+    ap.add_argument("--workers", type=int, default=16,
+                     help="number of quiet_filter.py subprocesses to run "
+                          "concurrently (default 16, kept below this box's "
+                          "24 cores to leave other jobs some breathing room)")
+    args = ap.parse_args()
+
+    with tempfile.TemporaryDirectory(prefix="shard_quiet_filter_") as tmpdir_str:
+        tmpdir = Path(tmpdir_str)
+        shard_ins = [tmpdir / f"shard_{i:02d}.in.jsonl" for i in range(args.workers)]
+        shard_outs = [tmpdir / f"shard_{i:02d}.out.jsonl" for i in range(args.workers)]
+
+        print(f"Sharding {args.infile} across {args.workers} workers ...",
+              file=sys.stderr)
+        shard_input(args.infile, shard_ins)
+
+        print(f"Launching {args.workers} quiet_filter.py workers ...",
+              file=sys.stderr)
+        procs = [
+            subprocess.Popen([sys.executable, str(_QUIET_FILTER),
+                               "--in", str(sin), "--out", str(sout)])
+            for sin, sout in zip(shard_ins, shard_outs)
+        ]
+        returncodes = [p.wait() for p in procs]
+        if any(rc != 0 for rc in returncodes):
+            sys.exit(f"one or more quiet_filter.py workers failed "
+                      f"(return codes: {returncodes})")
+
+        print(f"Combining {args.workers} shard outputs -> {args.outfile} ...",
+              file=sys.stderr)
+        n_rows = 0
+        with open(args.outfile, "w") as fout:
+            for sout in shard_outs:
+                with open(sout) as fin:
+                    for line in fin:
+                        fout.write(line)
+                        n_rows += 1
+
+    print(f"Wrote {args.outfile}: {n_rows} quiet positions "
+          f"(combined from {args.workers} shards)")
+
+
+if __name__ == "__main__":
+    main()