docs(research): SO-ISMCTS BC ceiling write-up from 2026-05-11 autonomous session

Summarizes the 13-cycle trap-exploration session: BC pretrain (heuristic
clone) is the self-play ceiling under our compute budget (1 GPU + 50
sims + 768x4 MLP). All variants (naive finetune, KL anchor, mirror
descent, mixed-opponent + opponent-aware search) either preserved BC
(~17-21/100 vs heuristic-cautious) or regressed to catastrophic
forgetting. The single largest improvement of the session — 4× win rate
on the same checkpoint — came from PUCT Q-value normalization at search
time, not from any learning change.

Records the mechanism (negative training signal from BC-vs-heuristic
games; search too shallow to find heuristic-beating moves), the
hypotheses we negated, and the dials left in code for future runs with
more compute.

configs/ismcts/default.yaml

@@ -1,8 +1,8 @@
 run:
   experiment_name: ismcts-default
-  max_iterations: 100
+  max_iterations: 500
   seed: 1
-  device: auto
+  device: cuda
 rules:
   n_colors: 5
   n_ranks: 9
@@ -17,15 +17,20 @@ encoding:
   slot_aware_playability: true
 network:
   kind: mlp
-  hidden_size: 512
+  hidden_size: 768
-  num_layers: 3
+  num_layers: 4
   activation: relu
 mcts:
   n_simulations: 50
-  c_puct: 1.5
+  c_puct: 5.0
   max_depth: 200
-  parallel_simulations: 8
+  parallel_simulations: 64
-  virtual_loss_value: 1.0
+  virtual_loss_value: 5.0
+  eval_n_simulations: 16
+  rollout_policy: heuristic_balanced
+  use_rollout_value: false
+  root_dirichlet_alpha: 0.3
+  root_dirichlet_epsilon: 0.4
 temperature:
   training: 1.0
   eval: 0.0
@@ -36,15 +41,17 @@ training:
   replay_capacity: 100000
   interleave_games: 8
   interleave_max_batch: 64
+  num_workers: 8
+  worker_device: cuda
 optimization:
   learning_rate: 0.0003
   grad_clip: 5.0
 checkpoint:
-  save_every: 10
+  save_every: 20
   save_latest: true
 evaluation:
-  eval_every: 10
+  eval_every: 5
-  games: 20
+  games: 5
   opponents: [random, discard-only, heuristic-cautious]
-  max_steps: 10000
+  max_steps: 500
-  num_workers: 1
+  num_workers: 8

docs/research/ismcts-bc-ceiling-2026-05-11.md

@@ -0,0 +1,169 @@
+# SO-ISMCTS BC Ceiling — 2026-05-11 Autonomous Session
+
+**Last verified:** 2026-05-11, commit `cba6cae` (branch `autonomous/trap-exploration`)
+
+## Short answer
+
+Under our current compute budget (1 GPU, 12 CPU cores, 50 MCTS sims/move,
+768x4 MLP), **behavior-cloning the heuristic-balanced bot is the ceiling**.
+Across 13 self-play training variants, no run cleared the BC baseline of
+21/100 wins vs `heuristic-cautious` in 100-game evaluation. Every variant
+either preserved BC (KL anchor, mirror-descent target) or regressed toward
+catastrophic forgetting (naive finetune, high-fraction mixed opponent).
+
+The single largest improvement of the session came from PUCT Q-value
+normalization at the *search* level, not from any learning change.
+
+## Headline numbers (vs heuristic-cautious, 100 games, n_sims = 16)
+
+| Run               | Setup                                          | W/100 | Notes |
+|-------------------|------------------------------------------------|------:|-------|
+| BC pretrain       | 5k heuristic-vs-heuristic games, 20 epochs CE+MSE |  21 | baseline |
+| C9 naive finetune | BC + plain self-play (no regularizer)          |   0 | catastrophic forgetting |
+| C10 KL β=1.0      | BC + self-play + KL(current ‖ BC)              | ~21 | preserved BC, no improvement |
+| C11 KL β=0.3      | weaker anchor                                  | ~21 | preserved BC, no improvement |
+| C12 mirror desc.  | target = softmax(α log π_mcts + (1−α) log π_BC) |  19 | preserved BC, no improvement |
+| C13 mixed=0.5     | 50 % games vs heuristic-balanced, opponent-aware MCTS, no KL | 0 | forgetting (worse than naive) |
+| C14 mixed=0.2     | mixed-opponent + opponent-aware + KL β=1.0     |  17 | preserved BC, no improvement |
+
+CIs (Wilson 95 %) overlap across all "preserved BC" rows; the 17–22 band
+is statistically indistinguishable from the BC baseline.
+
+## What actually moved the needle: PUCT Q normalization
+
+`mcts.pyx _select_action` previously used the raw score-unit Q:
+
+```
+score = q_eff + c_puct * prior * sqrt(N) / (1 + n)
+```
+
+With `value_scale = 100` (Lost Cities score units), a single backup could
+swing `q_eff` by ±100, while the exploration bonus is ~1–10. A noisy value
+at the root permanently buried low-prior actions before they could be
+explored.
+
+Fix (`b9fc569`): divide Q by `q_scale` (defaults to 100) before scoring:
+
+```
+score = q_eff / q_scale + c_puct * prior * sqrt(N) / (1 + n)
+```
+
+Replaying the exact same BC checkpoint with this fix took win rate vs
+heuristic-cautious from **5/100 → 21/100** — a 4× improvement from a
+~10-line search change, with no retraining. Worth holding onto as the
+load-bearing finding of the session.
+
+## Hypotheses we negated
+
+1. **Symmetric self-play eventually escapes the weak fixed point.**
+   Random-init + KL-free self-play ran for 100s of iterations across
+   C1–C8 without exceeding the noise floor (0–25 wins, all CIs overlap
+   each other and zero).
+
+2. **Mixed-opponent self-play (Codex top pick) breaks the weak
+   equilibrium.** With opponent-aware MCTS so the search distribution
+   reflects the real opponent (per Codex's "pitfall" warning), C13
+   regressed to 0/100. The training signal from vs-heuristic games is
+   structurally negative — BC cannot beat the heuristic, so every mixed
+   sample is a loss, and the gradient labels every BC action as bad.
+   C14 cut the fraction to 0.2 and added a strong KL anchor (β = 1.0),
+   which preserved BC but did not lift it.
+
+3. **Deeper search compensates for weak learning.** Increasing
+   `n_simulations` from 50 → 200 on the BC checkpoint *reduced* wins
+   vs `heuristic-balanced` from 28/64 → 14/64 in earlier probing.
+   Deeper search amplifies the network's preferences, including its
+   weaker ones, without supplying new information.
+
+4. **A different regularizer would let self-play improve on BC.**
+   KL anchor (β ∈ {0.3, 1.0}) and mirror-descent target mixing (α
+   annealed 0.3 → 0.8) both kept the network glued to BC. Neither
+   supplied a positive gradient to walk away from it.
+
+## Why BC is the ceiling — the mechanism
+
+Self-play seeded from a strong heuristic faces a structural trap:
+
+- BC has internalized the heuristic. Two BC copies playing each other
+  produce a near-symmetric outcome distribution; the visit counts at
+  most nodes give little policy-improvement signal beyond what BC
+  already encodes.
+- Against the real heuristic, BC loses systematically (the heuristic
+  beats its own clone in approx. 79 % of games at our scale). The
+  resulting training signal is uniformly negative; learning that
+  signal pushes the policy *away* from BC without pointing anywhere
+  productive.
+- With 50 MCTS sims/move on a 768x4 network, the search cannot
+  reliably *find* moves that beat the heuristic. So the only way out
+  of the trap — discovering a positive improvement direction —
+  is closed by the search-depth budget.
+
+The result is consistent with the standard SO-ISMCTS picture: π_weak
+(the symmetric weak fixed point) sits at roughly BC strength, π_Nash
+is unreachable at this compute, and every variant we tried collapses
+onto π_weak.
+
+## Things left as configurable dials (no behavior change at defaults)
+
+The `autonomous/trap-exploration` branch leaves the following in place
+for future runs with more compute:
+
+- `MctsConfig.q_scale` — PUCT Q normalization (defaults to 100, keep).
+- `MctsConfig.root_dirichlet_alpha / epsilon` — AlphaZero exploration noise.
+- `MctsConfig.opponent_aware_search` — when true, MCTS treats the
+  opponent seat as an external bot (skips tree expansion on opponent
+  turns, traverser-centered values).
+- `TrainingConfig.mixed_opponent_fraction` — 0 disables (pure self-play).
+- `TrainingConfig.mixed_opponent_bot` — bot name from
+  `coolrl_lost_cities.games.classic.bots.registry`.
+- `TrainingConfig.kl_anchor_ckpt` / `kl_anchor_beta` — frozen reference
+  network for `KL(current ‖ ref)` regularization.
+- `TrainingConfig.md_target_ref_ckpt` / `md_target_alpha_*` — mirror-
+  descent policy target with annealed mixing.
+- `lost-cities-ismcts pretrain` — heuristic behavior cloning subcommand.
+- `lost-cities-ismcts eval --ckpt … --n-sims N --games N --device cpu`
+  — standalone evaluator with Wilson CIs (`eval_checkpoint.py`).
+
+## What would be worth trying with more compute
+
+Not implemented here. These are the directions that the mechanism above
+*does not rule out*:
+
+- **Deeper search at training time** (n_sims ≫ 200, e.g. 800–1600).
+  Enough simulations should eventually surface a heuristic-beating
+  action somewhere in the search tree; that's a positive gradient.
+- **Population training with frozen snapshots.** Periodically snapshot
+  the trainer and route 10–20 % of self-play games against the snapshot
+  pool. Combined with opponent-aware search, this gives a stationary
+  diverse-opponent gradient without the all-negative-signal problem of
+  pure heuristic mixing.
+- **Value-weighted replay.** Prioritize high-error samples in the
+  buffer so the value head sees the cases where it disagrees with the
+  search rollout.
+- **Larger / better-shaped networks.** 768x4 MLP may simply lack the
+  capacity to represent the conjunctions Lost Cities needs (color ×
+  expedition × hand composition). Attention or factored heads could be
+  worth probing.
+
+## Code references
+
+- Search-side: `src/coolrl_lost_cities/games/classic/ismcts/mcts.pyx`
+  (`_select_action`, `prepare_simulation`, `_expand_with_prior`).
+- Python parity: `src/coolrl_lost_cities/games/classic/ismcts/mcts.py`.
+- Mixed-opponent / opponent-aware wiring:
+  `src/coolrl_lost_cities/games/classic/ismcts/interleaved_self_play.py`.
+- Regularization (KL anchor, mirror-descent) and metrics:
+  `src/coolrl_lost_cities/games/classic/ismcts/trainer.py`.
+- BC pretrain: `src/coolrl_lost_cities/games/classic/ismcts/pretrain.py`.
+- Eval CLI with Wilson CIs:
+  `src/coolrl_lost_cities/games/classic/ismcts/eval_checkpoint.py`.
+- BC checkpoint (load with `--resume-from`):
+  `runs/pretrain/heuristic_balanced_5kg_20ep.pt` (5 k games, 20 epochs).
+
+## Related memory
+
+- `opponent-policy-network-divergence.md` — the Deep CFR analogue:
+  using the live network as its own opponent breaks stationarity and
+  diverges. The SO-ISMCTS picture here is the same family of failure:
+  bootstrapping from oneself does not provide a positive learning
+  signal.

scripts/autonomous_cycle_eval.sh

@@ -0,0 +1,30 @@
+#!/bin/bash
+# Autonomous cycle eval helper.
+# Usage: ./autonomous_cycle_eval.sh <run-prefix> [extra eval args...]
+# Finds latest run matching prefix, runs eval --ckpt latest.pt with 30 games,
+# and reports: timeouts, natural-end wins per opponent.
+
+set -euo pipefail
+
+PREFIX="${1:-}"
+shift || true
+
+if [ -z "$PREFIX" ]; then
+    echo "usage: $0 <run-prefix> [extra eval args...]"
+    exit 1
+fi
+
+RUN=$(ls -td runs/*${PREFIX}* 2>/dev/null | head -1)
+if [ -z "$RUN" ]; then
+    echo "no run matching ${PREFIX}" >&2
+    exit 1
+fi
+
+CKPT="$RUN/latest.pt"
+if [ ! -f "$CKPT" ]; then
+    echo "no checkpoint at $CKPT" >&2
+    exit 1
+fi
+
+echo "=== eval $CKPT ==="
+uv run lost-cities-ismcts eval --ckpt "$CKPT" --games 30 --verbose "$@" 2>&1

src/coolrl_lost_cities/games/classic/ismcts/cli.py

@@ -79,6 +79,19 @@ def train_command(args: argparse.Namespace) -> None:
         device=config.run.device,
         tracker=tracker,
     )
+    if args.resume_from:
+        import torch
+
+        ckpt = torch.load(args.resume_from, map_location=trainer.device, weights_only=False)
+        trainer.network.load_state_dict(ckpt["network"])
+        if "optimizer" in ckpt:
+            trainer.optimizer.load_state_dict(ckpt["optimizer"])
+        print(
+            f"[resume] loaded network + optimizer from {args.resume_from} "
+            f"(prior iteration={ckpt.get('iteration', '?')}); "
+            f"new run starts at iteration 1 with current config",
+            flush=True,
+        )
     try:
         trainer.train()
     finally:
@@ -111,7 +124,31 @@ def main(argv: list[str] | None = None) -> None:
     train.add_argument("--wandb-job-type")
     train.add_argument("--wandb-tag", action="append", default=[])
     train.add_argument("--wandb-notes")
+    train.add_argument(
+        "--resume-from",
+        default=None,
+        help="Path to a .pt checkpoint to warm-start network + optimizer state.",
+    )
     train.set_defaults(func=train_command)
+
+    from .eval_checkpoint import add_eval_args, run_eval
+
+    eval_cmd = subparsers.add_parser(
+        "eval",
+        help="Evaluate a saved checkpoint vs heuristic bots in parallel.",
+    )
+    add_eval_args(eval_cmd)
+    eval_cmd.set_defaults(func=lambda a: run_eval(a))
+
+    from .pretrain import add_pretrain_args, run_pretrain
+
+    pretrain_cmd = subparsers.add_parser(
+        "pretrain",
+        help="Behavior-clone a heuristic bot into the AlphaZero network as a warm start.",
+    )
+    add_pretrain_args(pretrain_cmd)
+    pretrain_cmd.set_defaults(func=lambda a: run_pretrain(a))
+
     args = parser.parse_args(argv)
     args.func(args)
 

src/coolrl_lost_cities/games/classic/ismcts/config.py

@@ -30,6 +30,20 @@ class MctsConfig(StrictModel):
     virtual_loss_value: float = 1.0
     eval_with_mcts: bool = True
     eval_n_simulations: int = 0
+    root_dirichlet_alpha: float = 0.0
+    root_dirichlet_epsilon: float = 0.0
+    # Divisor applied to Q values inside PUCT to bring them onto roughly the
+    # same scale as the exploration bonus. With value_scale=100 score units,
+    # raw Q can swing ±100 while c_puct * prior * sqrt(N) is ~1-10, so a single
+    # bad backup permanently kills an action. Setting q_scale=100 normalizes Q
+    # to ~[-1, 1] (consistent with AlphaZero's convention).
+    q_scale: float = 100.0
+    # Opponent-aware search: when set, the search tree treats the opponent
+    # seat as a fixed external policy (heuristic bot) instead of expanding it
+    # with the network's priors/value. Used during mixed-opponent self-play
+    # so root visit distributions reflect the *actual* opponent the trainee
+    # faces. Bot name is taken from training.mixed_opponent_bot.
+    opponent_aware_search: bool = False
 
     @field_validator("n_simulations", "max_depth", "parallel_simulations")
     @classmethod
@@ -60,6 +74,37 @@ class TrainingConfig(StrictModel):
     interleave_max_batch: int = 64
     num_workers: int = 1
     worker_device: str = "cpu"
+    # Multiplier on the value-head MSE loss (already normalized by value_scale**2).
+    # Default 1.0 keeps current behavior; raising it (e.g. 50-100) makes the value
+    # head learn faster relative to policy loss. Useful when value_prediction_error
+    # is large but loss/value is tiny because of the normalization.
+    value_loss_weight: float = 1.0
+    # Optional KL anchor to a reference (e.g. behavior-cloned) policy. The
+    # reference network is loaded once at trainer start and frozen; on every
+    # gradient step we add `kl_anchor_beta * KL(current || reference)` to the
+    # loss. Anchors self-play training to the pretrained policy and prevents
+    # catastrophic forgetting / drift to weak self-play equilibria.
+    kl_anchor_ckpt: str | None = None
+    kl_anchor_beta: float = 0.0
+    # Mirror-descent target mixing for policy loss. Alternative to kl_anchor;
+    # blends MCTS visit distribution with the reference (BC) policy in log
+    # space, then trains the network to match. pi_target = softmax(
+    #   alpha * log(pi_mcts) + (1 - alpha) * log(pi_ref)
+    # ). Anneal alpha from low (rely on BC) to high (rely on MCTS) over
+    # training. Requires kl_anchor_ckpt to be set as the reference source.
+    md_target_ref_ckpt: str | None = None
+    md_target_alpha_start: float = 0.3
+    md_target_alpha_end: float = 0.8
+    md_target_alpha_iters: int = 500
+    # Mixed-opponent self-play: a fraction of games per iteration are played
+    # against a fixed external bot instead of the current network. Only the
+    # trainee's decisions are stored as policy targets; opponent moves are
+    # taken by `mixed_opponent_bot.act(state)`. Combined with
+    # mcts.opponent_aware_search, the MCTS tree models the opponent as that
+    # same bot so root-visit distributions reflect the real opponent.
+    # Set fraction=0 to disable (pure self-play).
+    mixed_opponent_fraction: float = 0.0
+    mixed_opponent_bot: str = "heuristic-balanced"
 
     @field_validator(
         "games_per_iter",

src/coolrl_lost_cities/games/classic/ismcts/eval_checkpoint.py

@@ -0,0 +1,333 @@
+"""Standalone parallel evaluation of an ISMCTS checkpoint vs heuristic bots.
+
+Reuses the same MCTS / bot stack as the training-loop eval, but does not
+interfere with a running training process. Useful for comparing a snapshot
+against opponents that are not in `evaluation.opponents` (e.g. heuristic-balanced,
+the rollout policy) and for running many more games than per-iter eval typically
+allows.
+
+Invoked via ``lost-cities-ismcts eval`` (see ``cli.py``).
+"""
+
+from __future__ import annotations
+
+import argparse
+import multiprocessing as mp
+import os
+import random
+import sys
+import time
+from concurrent.futures import ProcessPoolExecutor, as_completed
+from dataclasses import dataclass, field
+from pathlib import Path
+
+import torch
+
+
+@dataclass(frozen=True)
+class _WorkerJob:
+    ckpt_path: str
+    opponent: str
+    game_indices: tuple[int, ...]
+    seed: int
+    device: str
+    worker_index: int
+    verbose: bool
+    n_sims_override: int = 0  # 0 means use checkpoint's eval_n_simulations
+
+
+@dataclass
+class _GameResult:
+    game_index: int
+    policy_player: int
+    score_diff: float
+    turns: int
+    policy_turns: int
+    play_actions: int
+    timed_out: bool
+
+
+@dataclass
+class _WorkerResult:
+    worker_index: int
+    opponent: str
+    games: list[_GameResult] = field(default_factory=list)
+    elapsed: float = 0.0
+
+
+def _run_games(job: _WorkerJob) -> _WorkerResult:
+    # Inside-worker imports + thread-cap to avoid CPU oversubscription
+    from coolrl_lost_cities.games.classic.bots.registry import build_bot
+    from coolrl_lost_cities.games.classic.deep_cfr.encoding import input_dim
+    from coolrl_lost_cities.games.classic.game import GameState, LostCitiesConfig
+    from coolrl_lost_cities.games.classic.ismcts.config import IsMctsConfig
+    from coolrl_lost_cities.games.classic.ismcts.mcts import IsMctsSearcher
+    from coolrl_lost_cities.games.classic.ismcts.network import AlphaZeroNet
+
+    os.environ.setdefault("OMP_NUM_THREADS", "1")
+    os.environ.setdefault("MKL_NUM_THREADS", "1")
+    torch.set_num_threads(1)
+
+    started = time.perf_counter()
+    print(
+        f"  [worker {job.worker_index}] start ({len(job.game_indices)} games "
+        f"vs {job.opponent}, device={job.device})",
+        flush=True,
+    )
+
+    ckpt = torch.load(job.ckpt_path, map_location="cpu", weights_only=False)
+    cfg = IsMctsConfig.model_validate(ckpt["config"])
+    game_config = LostCitiesConfig(**ckpt["game_config"])
+    probe = GameState.new_game(game_config, seed=cfg.run.seed)
+    dim = input_dim(probe, cfg.encoding)
+    device = torch.device(job.device)
+    net = AlphaZeroNet.from_config(dim, probe.action_size, cfg).to(device)
+    net.load_state_dict(ckpt["network"])
+    net.eval()
+
+    eval_mcts_cfg = cfg.mcts.model_copy()
+    if job.n_sims_override > 0:
+        eval_mcts_cfg = eval_mcts_cfg.model_copy(update={"n_simulations": job.n_sims_override})
+    elif cfg.mcts.eval_n_simulations > 0:
+        eval_mcts_cfg = eval_mcts_cfg.model_copy(
+            update={"n_simulations": cfg.mcts.eval_n_simulations}
+        )
+
+    rng = random.Random(job.seed + job.worker_index * 7919)
+    result = _WorkerResult(worker_index=job.worker_index, opponent=job.opponent)
+
+    max_steps = 500
+    for game_index in job.game_indices:
+        game_started = time.perf_counter()
+        policy_player = game_index % 2
+        opps = [
+            build_bot(job.opponent, seed=job.seed + game_index),
+            build_bot(job.opponent, seed=job.seed + game_index + 1),
+        ]
+        state = GameState.new_game(game_config, seed=job.seed + game_index)
+        turns = 0
+        policy_turns = 0
+        play_actions = 0
+        timed_out = False
+        while True:
+            if state.terminal:
+                break
+            if turns >= max_steps:
+                timed_out = True
+                break
+            current = int(state.current_player)
+            if current == policy_player:
+                searcher = IsMctsSearcher(
+                    net,
+                    eval_mcts_cfg,
+                    device=device,
+                    encoding=cfg.encoding,
+                    rng=random.Random(rng.randrange(2**31)),
+                )
+                visits = searcher.search(state, current)
+                unified = (
+                    max(visits, key=visits.get) if visits else state.unified_legal_actions()[0]
+                )
+                if state.phase == "card":
+                    policy_turns += 1
+                    if unified % 2 == 0:
+                        play_actions += 1
+                state.apply_unified_action(unified)
+            else:
+                state.apply_action(opps[current].act(state))
+            turns += 1
+
+        diff = float(state.score_diff(policy_player))
+        gr = _GameResult(
+            game_index=game_index,
+            policy_player=policy_player,
+            score_diff=diff,
+            turns=turns,
+            policy_turns=policy_turns,
+            play_actions=play_actions,
+            timed_out=timed_out,
+        )
+        result.games.append(gr)
+        if job.verbose:
+            elapsed_g = time.perf_counter() - game_started
+            pa = play_actions / policy_turns if policy_turns else 0.0
+            print(
+                f"  [worker {job.worker_index}] game {game_index:3d} "
+                f"as P{policy_player} | turns={turns:3d} "
+                f"score={diff:+6.1f} PA={pa:.2f}"
+                f"{' TIMEOUT' if timed_out else ''} "
+                f"({elapsed_g:.1f}s)",
+                flush=True,
+            )
+
+    result.elapsed = time.perf_counter() - started
+    won = sum(1 for g in result.games if g.score_diff > 0)
+    print(
+        f"  [worker {job.worker_index}] done {len(result.games)} games "
+        f"vs {job.opponent} in {result.elapsed:.1f}s "
+        f"(W={won}/{len(result.games)})",
+        flush=True,
+    )
+    return result
+
+
+def _split_games(n_games: int, n_workers: int) -> list[tuple[int, ...]]:
+    n_workers = max(1, min(n_workers, n_games))
+    base = n_games // n_workers
+    rem = n_games % n_workers
+    out: list[tuple[int, ...]] = []
+    cursor = 0
+    for i in range(n_workers):
+        count = base + (1 if i < rem else 0)
+        out.append(tuple(range(cursor, cursor + count)))
+        cursor += count
+    return out
+
+
+def _summarize(games: list[_GameResult]) -> dict[str, float]:
+    n = len(games)
+    if n == 0:
+        return {}
+    wins = sum(1 for g in games if g.score_diff > 0)
+    losses = sum(1 for g in games if g.score_diff < 0)
+    draws = sum(1 for g in games if g.score_diff == 0)
+    timeouts = sum(1 for g in games if g.timed_out)
+    score_diffs = [g.score_diff for g in games]
+    avg = sum(score_diffs) / n
+    var = sum((d - avg) ** 2 for d in score_diffs) / n if n > 1 else 0.0
+    std = var**0.5
+    total_policy_turns = sum(g.policy_turns for g in games)
+    total_play_actions = sum(g.play_actions for g in games)
+    pa = total_play_actions / total_policy_turns if total_policy_turns else 0.0
+    avg_turns = sum(g.turns for g in games) / n
+    z = 1.96
+    # Wilson CI on win rate (half-width only; report as wr ± half)
+    wr = wins / n
+    denom = 1 + z * z / n
+    half = z / denom * ((wr * (1 - wr) / n + z * z / (4 * n * n)) ** 0.5)
+    score_ci = z * std / (n**0.5) if n > 1 else 0.0
+    return {
+        "games": n,
+        "wins": wins,
+        "losses": losses,
+        "draws": draws,
+        "timeouts": timeouts,
+        "win_rate": wr,
+        "win_rate_ci_half": half,
+        "avg_score_diff": avg,
+        "score_std": std,
+        "score_ci_half": score_ci,
+        "play_action_rate": pa,
+        "avg_turns": avg_turns,
+    }
+
+
+def _default_ckpt() -> Path:
+    """Find latest 'ismcts-overnight*' run's latest.pt, or fallback to newest run."""
+    candidates = sorted(Path("runs").glob("*ismcts-overnight*"))
+    if candidates:
+        return candidates[-1] / "latest.pt"
+    candidates = sorted(Path("runs").iterdir())
+    if not candidates:
+        raise SystemExit("no runs/ directory entries")
+    return candidates[-1] / "latest.pt"
+
+
+def add_eval_args(parser: argparse.ArgumentParser) -> None:
+    parser.add_argument(
+        "--ckpt",
+        default=None,
+        help="Path to checkpoint .pt. Default: latest overnight run latest.pt.",
+    )
+    parser.add_argument(
+        "--opponents",
+        nargs="+",
+        default=["heuristic-balanced", "heuristic-aggressive", "heuristic-cautious"],
+        help="Opponent bot names from registry.",
+    )
+    parser.add_argument("--games", type=int, default=50, help="Games per opponent (default: 50).")
+    parser.add_argument(
+        "--device",
+        choices=("cpu", "cuda"),
+        default="cpu",
+        help="Worker device (default: cpu; cuda may compete with running training).",
+    )
+    parser.add_argument(
+        "--num-workers", type=int, default=8, help="Parallel worker count (default: 8)."
+    )
+    parser.add_argument("--seed", type=int, default=99999)
+    parser.add_argument(
+        "--verbose",
+        action="store_true",
+        help="Print per-game result lines (turns/score/PA) in addition to per-worker summaries.",
+    )
+    parser.add_argument(
+        "--n-sims",
+        type=int,
+        default=0,
+        help="Override n_simulations at eval time (0 = use checkpoint's eval_n_simulations).",
+    )
+
+
+def run_eval(args: argparse.Namespace) -> None:
+    ckpt_path = Path(args.ckpt) if args.ckpt else _default_ckpt()
+    if not ckpt_path.exists():
+        print(f"checkpoint not found: {ckpt_path}", file=sys.stderr)
+        raise SystemExit(1)
+
+    ckpt = torch.load(ckpt_path, map_location="cpu", weights_only=False)
+    iteration = ckpt.get("iteration", "?")
+    print(f"checkpoint : {ckpt_path}")
+    print(f"iteration  : {iteration}")
+    print(f"device     : {args.device}")
+    print(f"workers    : {args.num_workers}")
+    print(f"games/opp  : {args.games}")
+    print(f"opponents  : {args.opponents}")
+    print(f"seed       : {args.seed}")
+    if args.verbose:
+        print("verbose    : True (per-game logging)")
+    print()
+
+    ctx = mp.get_context("spawn")
+    started_all = time.perf_counter()
+
+    for opponent in args.opponents:
+        opp_started = time.perf_counter()
+        slices = _split_games(args.games, args.num_workers)
+        jobs = [
+            _WorkerJob(
+                ckpt_path=str(ckpt_path),
+                opponent=opponent,
+                game_indices=tuple(slices[i]),
+                seed=args.seed,
+                device=args.device,
+                worker_index=i,
+                verbose=args.verbose,
+                n_sims_override=int(args.n_sims),
+            )
+            for i in range(len(slices))
+        ]
+        all_games: list[_GameResult] = []
+        with ProcessPoolExecutor(max_workers=len(jobs), mp_context=ctx) as ex:
+            futures = [ex.submit(_run_games, j) for j in jobs]
+            for f in as_completed(futures):
+                res = f.result()
+                all_games.extend(res.games)
+
+        elapsed = time.perf_counter() - opp_started
+        summary = _summarize(all_games)
+        n = int(summary["games"])
+        print()
+        print(
+            f"vs {opponent:22s} | W={summary['wins']}/{n} "
+            f"({summary['win_rate']:.2f} ± {summary['win_rate_ci_half']:.2f}) "
+            f"| S={summary['avg_score_diff']:+6.1f} ± {summary['score_ci_half']:5.1f} "
+            f"(σ={summary['score_std']:.1f}) | PA={summary['play_action_rate']:.2f} "
+            f"| turns={summary['avg_turns']:.0f} "
+            f"| timeouts={summary['timeouts']} "
+            f"| elapsed={elapsed:.1f}s"
+        )
+        print()
+
+    total = time.perf_counter() - started_all
+    print(f"total elapsed: {total:.1f}s")

src/coolrl_lost_cities/games/classic/ismcts/interleaved_self_play.py

@@ -6,6 +6,7 @@ from dataclasses import dataclass, field
 import numpy as np
 import torch
 
+from coolrl_lost_cities.games.classic.bots.registry import build_bot
 from coolrl_lost_cities.games.classic.deep_cfr.encoding import encode_info_state
 from coolrl_lost_cities.games.classic.game import GameState, LostCitiesConfig
 
@@ -34,6 +35,11 @@ class _GameContext:
     game_index: int
     decisions: list[_PendingDecision] = field(default_factory=list)
     steps: int = 0
+    # Mixed-opponent setup: when traverser_seat is not None, only that seat
+    # uses MCTS+network; the other seat is played by `opponent_bot`. None for
+    # pure self-play games (both seats use MCTS).
+    traverser_seat: int | None = None
+    opponent_bot: object | None = None
 
 
 @dataclass
@@ -61,15 +67,28 @@ def play_self_play_iteration(
     active: list[_GameContext] = []
     started = 0
     target_games = training_config.games_per_iter
+    mixed_fraction = float(training_config.mixed_opponent_fraction)
 
     def fill_active() -> None:
         nonlocal started
         while len(active) < training_config.interleave_games and started < target_games:
+            traverser_seat: int | None = None
+            opponent_bot = None
+            if mixed_fraction > 0.0 and rng.random() < mixed_fraction:
+                # Alternate trainee seat so MCTS sees both first- and
+                # second-player perspectives equally.
+                traverser_seat = started % 2
+                opponent_bot = build_bot(
+                    training_config.mixed_opponent_bot,
+                    seed=rng.randrange(2**31),
+                )
             active.append(
                 _GameContext(
                     state=GameState.new_game(game_config, seed=rng.randrange(2**31)),
                     rng=random.Random(rng.randrange(2**31)),
                     game_index=started,
+                    traverser_seat=traverser_seat,
+                    opponent_bot=opponent_bot,
                 )
             )
             started += 1
@@ -83,6 +102,19 @@ def play_self_play_iteration(
                 completed.append(_finalize_context(context))
                 continue
             player = int(context.state.current_player)
+            # Mixed-opponent: if it's the opponent's turn in a mixed game,
+            # let the heuristic bot move directly (no MCTS, no sample).
+            if (
+                context.traverser_seat is not None
+                and context.opponent_bot is not None
+                and player != context.traverser_seat
+            ):
+                phase_action = context.opponent_bot.act(context.state)
+                unified = context.state.to_unified_action(phase_action)
+                context.state.apply_unified_action(unified)
+                context.steps += 1
+                still_active.append(context)
+                continue
             searcher = IsMctsSearcher(
                 network,
                 mcts_config,
@@ -90,6 +122,18 @@ def play_self_play_iteration(
                 encoding=encoding,
                 rng=random.Random(context.rng.randrange(2**31)),
             )
+            # Pass opponent bot into the searcher for opponent-aware
+            # determinization (search models the real opponent the trainee
+            # faces, not a self-play mirror).
+            if (
+                mcts_config.opponent_aware_search
+                and context.traverser_seat is not None
+                and context.opponent_bot is not None
+            ):
+                searcher.set_opponent_bot(
+                    context.opponent_bot,
+                    traverser_seat=context.traverser_seat,
+                )
             jobs.append(
                 _SearchJob(
                     context=context,
@@ -181,6 +225,7 @@ def _evaluate_global_batch(
                 item.legal_actions,
                 priors_by_id[id(item)],
                 values_by_id[id(item)],
+                not item.path,  # is_root for Dirichlet noise
             )
         job.searcher._backup(item.path, value, item.leaf_player)
 
@@ -222,7 +267,15 @@ def _finish_decision(
 
 
 def _finalize_context(context: _GameContext) -> list[ReplaySample]:
-    final_diff0 = float(context.state.score_diff(0))
+    # If the game did not terminate naturally (hit max_steps), the score
+    # reflects an incomplete game — typically a "stall" outcome where both
+    # sides have under-developed expeditions. Treating that as a real win
+    # for either player creates a degenerate stall-and-pray learning
+    # signal. Zero it out so the trajectory is neutral.
+    if not context.state.terminal:
+        final_diff0 = 0.0
+    else:
+        final_diff0 = float(context.state.score_diff(0))
     samples: list[ReplaySample] = []
     for decision in context.decisions:
         value = final_diff0 if decision.player == 0 else -final_diff0

src/coolrl_lost_cities/games/classic/ismcts/mcts.py

@@ -89,6 +89,13 @@ class IsMctsSearcher:
         self._rollout_bot = (
             HeuristicBot() if config.rollout_policy == "heuristic_balanced" else None
         )
+        # Opponent-aware search: see mcts.pyx for the rationale.
+        self._opponent_bot: object | None = None
+        self._traverser_seat: int = -1
+
+    def set_opponent_bot(self, bot: object, *, traverser_seat: int) -> None:
+        self._opponent_bot = bot
+        self._traverser_seat = int(traverser_seat)
 
     def search(
         self,
@@ -133,34 +140,45 @@ class IsMctsSearcher:
         # Cache the info-set key for the current node so we don't recompute it
         # after applying an action (the child's key becomes the next iter's key).
         cached_key: bytes | None = None
+        opponent_aware = self._opponent_bot is not None
+        trav_seat = self._traverser_seat
         while True:
             player = int(state.current_player)
             if state.terminal or depth >= self.config.max_depth:
+                leaf_seat = trav_seat if opponent_aware else player
                 return PendingSimulation(
                     path=path,
                     leaf_state=state,
                     leaf_node=None,
-                    leaf_player=player,
+                    leaf_player=leaf_seat,
                     info_state=None,
                     legal_mask=None,
                     legal_actions=[],
-                    terminal_value=float(state.score_diff(player)),
+                    terminal_value=float(state.score_diff(leaf_seat)),
                 )
+            if opponent_aware and player != trav_seat:
+                phase_action = self._opponent_bot.act(state)
+                unified = state.to_unified_action(phase_action)
+                state.apply_unified_action(unified)
+                cached_key = None
+                depth += 1
+                continue
             key = cached_key if cached_key is not None else canonical_info_set_key(state, player)
             node = self.tree.get_or_create(key, player=player, terminal=state.terminal)
             if not node.is_expanded():
                 legal_actions = state.unified_legal_actions()
                 if not legal_actions:
                     node.terminal = True
+                    leaf_seat = trav_seat if opponent_aware else player
                     return PendingSimulation(
                         path=path,
                         leaf_state=state,
                         leaf_node=node,
-                        leaf_player=player,
+                        leaf_player=leaf_seat,
                         info_state=None,
                         legal_mask=None,
                         legal_actions=[],
-                        terminal_value=float(state.score_diff(player)),
+                        terminal_value=float(state.score_diff(leaf_seat)),
                     )
                 return PendingSimulation(
                     path=path,
@@ -263,6 +281,7 @@ class IsMctsSearcher:
         sqrt_total = math.sqrt(max(1, total_visits))
         best_score = -float("inf")
         best_action = legal_actions[0]
+        q_scale = float(getattr(self.config, "q_scale", 100.0)) or 1.0
         for action in legal_actions:
             n = node.visits.get(action, 0)
             virtual = node.virtual_visits.get(action, 0)
@@ -274,7 +293,8 @@ class IsMctsSearcher:
                 q_eff = (
                     node.value_sum.get(action, 0.0) - virtual * self.config.virtual_loss_value
                 ) / n_eff
-            score = q_eff + self.config.c_puct * prior * sqrt_total / (1 + n_eff)
+            # Normalize Q to match exploration-bonus scale; see mcts.pyx for details.
+            score = q_eff / q_scale + self.config.c_puct * prior * sqrt_total / (1 + n_eff)
             if score > best_score:
                 best_score = score
                 best_action = action

src/coolrl_lost_cities/games/classic/ismcts/mcts.pyx

@@ -296,6 +296,12 @@ cdef class IsMctsSearcher:
     cdef public MctsTree tree
     cdef HeuristicBot _rollout_bot
     cdef int action_size
+    # Opponent-aware search: when set, the search treats one seat as a fixed
+    # external policy (heuristic bot). Opponent moves are applied directly
+    # without entering the tree, and all values are taken from the
+    # traverser's perspective. None for standard symmetric self-play search.
+    cdef public object _opponent_bot
+    cdef public int _traverser_seat
 
     def __init__(
         self,
@@ -318,6 +324,12 @@ cdef class IsMctsSearcher:
         self._rollout_bot = (
             <HeuristicBot>PyHeuristicBot() if config.rollout_policy == "heuristic_balanced" else None
         )
+        self._opponent_bot = None
+        self._traverser_seat = -1
+
+    def set_opponent_bot(self, object bot, *, int traverser_seat):
+        self._opponent_bot = bot
+        self._traverser_seat = traverser_seat
 
     cdef inline int _from_unified_action_c(self, GameState state, int action_id) noexcept:
         cdef int card_action_size = 2 * state.hand_size
@@ -338,12 +350,16 @@ cdef class IsMctsSearcher:
         )
         cdef int sims = int(n_sims or self.config.n_simulations)
         cdef int completed = 0
+        cdef int batch_size
         cdef list pending
         cdef list legal
         cdef int action
         cdef dict result
         while completed < sims:
-            pending = self.prepare_simulation_batch(state, traverser, 1)
+            batch_size = min(int(self.config.parallel_simulations), sims - completed)
+            if batch_size <= 0:
+                batch_size = 1
+            pending = self.prepare_simulation_batch(state, traverser, batch_size)
             if not pending:
                 break
             self.evaluate_and_backup(pending)
@@ -386,19 +402,41 @@ cdef class IsMctsSearcher:
         cdef int actions[MAX_ACTIONS]
         cdef int action_count
         cdef int i
+        cdef bint opponent_aware = self._opponent_bot is not None
+        cdef int leaf_player_seat
+        cdef int trav_seat = self._traverser_seat
+        cdef object phase_action
+        cdef int unified_action
         while True:
             player = state.current_player
             if state.terminal or depth >= int(self.config.max_depth):
+                if opponent_aware:
+                    leaf_player_seat = trav_seat
+                else:
+                    leaf_player_seat = player
                 return PendingSimulation(
                     path=path,
                     leaf_state=state,
                     leaf_node=None,
-                    leaf_player=player,
+                    leaf_player=leaf_player_seat,
                     info_state=None,
                     legal_mask=None,
                     legal_actions=[],
-                    terminal_value=float(state.total_scores[player] - state.total_scores[1 - player]),
+                    terminal_value=float(
+                        state.total_scores[leaf_player_seat]
+                        - state.total_scores[1 - leaf_player_seat]
+                    ),
                 )
+            # Opponent-aware: if it's the opponent's turn, let the heuristic
+            # bot move directly instead of expanding the tree.
+            if opponent_aware and player != trav_seat:
+                phase_action = self._opponent_bot.act(state)
+                unified_action = state.to_unified_action(phase_action)
+                local_action = self._from_unified_action_c(state, unified_action)
+                state._push_action_c(local_action)
+                cached_key = None
+                depth += 1
+                continue
             if cached_key is None:
                 key = canonical_info_set_key(state, player)
             else:
@@ -409,15 +447,22 @@ cdef class IsMctsSearcher:
                 legal_actions = [actions[i] for i in range(action_count)]
                 if not legal_actions:
                     node.terminal = True
+                    if opponent_aware:
+                        leaf_player_seat = trav_seat
+                    else:
+                        leaf_player_seat = player
                     return PendingSimulation(
                         path=path,
                         leaf_state=state,
                         leaf_node=node,
-                        leaf_player=player,
+                        leaf_player=leaf_player_seat,
                         info_state=None,
                         legal_mask=None,
                         legal_actions=[],
-                        terminal_value=float(state.total_scores[player] - state.total_scores[1 - player]),
+                        terminal_value=float(
+                            state.total_scores[leaf_player_seat]
+                            - state.total_scores[1 - leaf_player_seat]
+                        ),
                     )
                 return PendingSimulation(
                     path=path,
@@ -489,6 +534,7 @@ cdef class IsMctsSearcher:
                     item.legal_actions,
                     priors_by_id[id(item)],
                     values_by_id[id(item)],
+                    not item.path,  # is_root: empty path means leaf == root
                 )
             self._backup(item.path, value, item.leaf_player)
 
@@ -500,14 +546,32 @@ cdef class IsMctsSearcher:
         list legal_actions,
         object probs,
         double network_value,
+        bint is_root=False,
     ):
         cdef int action
+        cdef int i
+        cdef int n_legal
+        cdef double alpha
+        cdef double epsilon
         cdef object rollout_value
+        cdef object noise
+        cdef object np_rng
         legal_actions = self._unified_legal_actions_list_c(state)
         if not legal_actions:
             node.terminal = True
             return float(state.total_scores[player] - state.total_scores[1 - player])
         node.expanded = True
+        # Dirichlet noise at root (AlphaZero pattern: force exploration of low-prior actions)
+        alpha = float(self.config.root_dirichlet_alpha)
+        epsilon = float(self.config.root_dirichlet_epsilon)
+        if is_root and alpha > 0.0 and epsilon > 0.0:
+            n_legal = len(legal_actions)
+            # Use numpy with seed from self.rng for reproducibility under fixed seeds
+            np_rng = np.random.default_rng(self.rng.randrange(2**31))
+            noise = np_rng.dirichlet([alpha] * n_legal)
+            for i in range(n_legal):
+                action = legal_actions[i]
+                probs[action] = (1.0 - epsilon) * float(probs[action]) + epsilon * float(noise[i])
         for action in legal_actions:
             (<_ArrayMap>node.priors).set_float(action, float(probs[action]))
             if not (<_ArrayMap>node.visits).has(action):
@@ -530,13 +594,17 @@ cdef class IsMctsSearcher:
         cdef double sqrt_total
         cdef double prior
         cdef double q_eff
+        cdef double q_normalized
         cdef double score
         cdef double best_score = -float("inf")
         cdef int best_action = int(legal_actions[0])
+        cdef double q_scale = float(getattr(self.config, "q_scale", 100.0))
         cdef _ArrayMap visits = <_ArrayMap>node.visits
         cdef _ArrayMap virtual_visits = <_ArrayMap>node.virtual_visits
         cdef _ArrayMap priors = <_ArrayMap>node.priors
         cdef _ArrayMap value_sum = <_ArrayMap>node.value_sum
+        if q_scale <= 0.0:
+            q_scale = 1.0
         for action in legal_actions:
             total_visits += visits.get_int(action, 0) + virtual_visits.get_int(action, 0)
         sqrt_total = math.sqrt(max(1, total_visits))
@@ -552,7 +620,12 @@ cdef class IsMctsSearcher:
                     value_sum.get_float(action, 0.0)
                     - virtual * float(self.config.virtual_loss_value)
                 ) / n_eff
-            score = q_eff + float(self.config.c_puct) * prior * sqrt_total / (1 + n_eff)
+            # Normalize Q to roughly [-1, 1] so the exploration bonus
+            # (c_puct * prior * sqrt(N) / (1+n)) competes on the right scale.
+            # Without this, raw score-units Q (±100) dominates and a single
+            # noisy backup kills exploration of low-prior actions.
+            q_normalized = q_eff / q_scale
+            score = q_normalized + float(self.config.c_puct) * prior * sqrt_total / (1 + n_eff)
             if score > best_score:
                 best_score = score
                 best_action = action

src/coolrl_lost_cities/games/classic/ismcts/pretrain.py

@@ -0,0 +1,279 @@
+"""Behavior cloning warm-start for the SO-ISMCTS network.
+
+Generates games from a fixed heuristic policy vs itself, then trains the
+AlphaZero-style network's policy + value heads in supervised fashion:
+
+- policy loss: cross-entropy between network logits and the heuristic's chosen
+  action (one-hot target, masked to legal actions).
+- value loss: MSE on the final game score diff from each decision-maker's
+  perspective, normalized by value_scale (same convention as the trainer).
+
+The resulting checkpoint can be passed to `lost-cities-ismcts train
+--resume-from` to start self-play with a heuristic-level prior instead of a
+random init. This addresses the self-play "weak-equilibrium" problem: starting
+from random, MCTS visit distributions converge to a mutually mediocre policy
+that has near-zero win rate against the heuristic. Warm-starting at heuristic
+level gives self-play a meaningful baseline to improve from.
+
+Invoked via ``lost-cities-ismcts pretrain``.
+"""
+
+from __future__ import annotations
+
+import argparse
+import random
+import time
+from pathlib import Path
+
+import numpy as np
+import torch
+import torch.nn.functional as F
+from torch import nn
+
+from coolrl_lost_cities.games.classic.bots.registry import build_bot
+from coolrl_lost_cities.games.classic.deep_cfr.encoding import encode_info_state, input_dim
+from coolrl_lost_cities.games.classic.game import GameState, LostCitiesConfig
+
+from .config import IsMctsConfig, load_config
+from .network import AlphaZeroNet
+
+
+def _collect_samples(
+    game_config: LostCitiesConfig,
+    encoding,
+    n_games: int,
+    bot_name: str = "heuristic-balanced",
+    seed: int = 0,
+    max_turns: int = 500,
+) -> list[tuple[np.ndarray, np.ndarray, int, float]]:
+    """Roll out n_games of bot vs bot, returning per-decision samples.
+
+    Each sample: (info_state, legal_mask, action_idx, value).
+    value is the player-perspective score diff at game end.
+    """
+    samples: list[tuple[np.ndarray, np.ndarray, int, float]] = []
+    for game_idx in range(n_games):
+        bots = [
+            build_bot(bot_name, seed=seed + game_idx * 2),
+            build_bot(bot_name, seed=seed + game_idx * 2 + 1),
+        ]
+        state = GameState.new_game(game_config, seed=seed + game_idx)
+        decisions: list[tuple[np.ndarray, np.ndarray, int, int]] = []
+        turns = 0
+        while not state.terminal and turns < max_turns:
+            player = int(state.current_player)
+            info_state = encode_info_state(state, player, encoding)
+            legal_mask = np.asarray(state.unified_legal_mask(), dtype=bool)
+            phase_action = bots[player].act(state)
+            unified = state.to_unified_action(phase_action)
+            decisions.append((info_state, legal_mask, int(unified), player))
+            state.apply_unified_action(unified)
+            turns += 1
+        final_diff0 = float(state.score_diff(0))
+        for info, mask, act, player in decisions:
+            value = final_diff0 if player == 0 else -final_diff0
+            samples.append((info, mask, act, value))
+    return samples
+
+
+def _train_supervised(
+    network: AlphaZeroNet,
+    samples: list,
+    device: torch.device,
+    *,
+    epochs: int,
+    batch_size: int,
+    lr: float,
+    weight_decay: float,
+    grad_clip: float,
+    value_loss_weight: float,
+) -> None:
+    network.train()
+    optimizer = torch.optim.AdamW(network.parameters(), lr=lr, weight_decay=max(weight_decay, 1e-4))
+    rng = random.Random(0)
+    indices = list(range(len(samples)))
+    v_scale = float(network.value_scale)
+
+    for epoch in range(1, epochs + 1):
+        rng.shuffle(indices)
+        n_batches = (len(indices) + batch_size - 1) // batch_size
+        epoch_pl = 0.0
+        epoch_vl = 0.0
+        epoch_acc = 0.0
+        epoch_n = 0
+        for b in range(n_batches):
+            batch_idx = indices[b * batch_size : (b + 1) * batch_size]
+            infos = np.stack([samples[i][0] for i in batch_idx])
+            masks = np.stack([samples[i][1] for i in batch_idx])
+            actions = np.array([samples[i][2] for i in batch_idx], dtype=np.int64)
+            values = np.array([samples[i][3] for i in batch_idx], dtype=np.float32)
+
+            info_t = torch.as_tensor(infos, dtype=torch.float32, device=device)
+            mask_t = torch.as_tensor(masks, dtype=torch.bool, device=device)
+            action_t = torch.as_tensor(actions, device=device)
+            value_t = torch.as_tensor(values, device=device)
+
+            logits, value_pred = network(info_t, mask_t)
+            policy_loss = F.cross_entropy(logits, action_t)
+            value_loss = F.mse_loss(value_pred / v_scale, value_t / v_scale)
+            loss = policy_loss + value_loss_weight * value_loss
+
+            optimizer.zero_grad(set_to_none=True)
+            loss.backward()
+            if grad_clip > 0:
+                nn.utils.clip_grad_norm_(network.parameters(), grad_clip)
+            optimizer.step()
+
+            with torch.no_grad():
+                preds = logits.argmax(dim=-1)
+                acc = (preds == action_t).float().mean().item()
+                epoch_pl += float(policy_loss.item()) * len(batch_idx)
+                epoch_vl += float(value_loss.item()) * len(batch_idx)
+                epoch_acc += acc * len(batch_idx)
+                epoch_n += len(batch_idx)
+
+        print(
+            f"  epoch {epoch:3d}: policy_loss={epoch_pl / epoch_n:.4f} "
+            f"value_loss={epoch_vl / epoch_n:.4f} "
+            f"top1_match={epoch_acc / epoch_n:.3f}",
+            flush=True,
+        )
+
+    network.eval()
+    return optimizer
+
+
+def add_pretrain_args(parser: argparse.ArgumentParser) -> None:
+    parser.add_argument(
+        "--config", default=None, help="ISMCTS config YAML (controls network shape + rules)."
+    )
+    parser.add_argument(
+        "--set",
+        action="append",
+        default=[],
+        dest="config_overrides",
+        metavar="PATH=VALUE",
+    )
+    parser.add_argument(
+        "--bot",
+        default="heuristic-balanced",
+        help="Bot to clone (default: heuristic-balanced).",
+    )
+    parser.add_argument(
+        "--games", type=int, default=2000, help="Number of bot-vs-bot games to roll out."
+    )
+    parser.add_argument("--epochs", type=int, default=10, help="Supervised training epochs.")
+    parser.add_argument("--batch-size", type=int, default=256)
+    parser.add_argument("--lr", type=float, default=3.0e-4)
+    parser.add_argument("--weight-decay", type=float, default=1.0e-4)
+    parser.add_argument("--grad-clip", type=float, default=5.0)
+    parser.add_argument(
+        "--value-loss-weight",
+        type=float,
+        default=50.0,
+        help="Multiplier on value MSE (raw_MSE / value_scale^2). Default 50 to "
+        "make value loss magnitude comparable to policy CE.",
+    )
+    parser.add_argument(
+        "--out",
+        default="runs/pretrain/heuristic_clone.pt",
+        help="Output checkpoint path (compatible with --resume-from).",
+    )
+    parser.add_argument("--device", default="cuda")
+    parser.add_argument("--seed", type=int, default=12345)
+
+
+def _apply_config_overrides(config: IsMctsConfig, assignments: list[str]) -> IsMctsConfig:
+    import yaml
+
+    def deep_update(base: dict, patch: dict) -> None:
+        for k, v in patch.items():
+            if isinstance(v, dict) and isinstance(base.get(k), dict):
+                deep_update(base[k], v)
+            else:
+                base[k] = v
+
+    overrides: dict = {}
+    for assignment in assignments:
+        if "=" not in assignment:
+            raise ValueError(f"override must be PATH=VALUE: {assignment}")
+        path, raw_value = assignment.split("=", 1)
+        value = yaml.safe_load(raw_value)
+        cursor = overrides
+        keys = path.split(".")
+        for k in keys[:-1]:
+            cursor = cursor.setdefault(k, {})
+        cursor[keys[-1]] = value
+    data = config.model_dump(mode="python")
+    deep_update(data, overrides)
+    return IsMctsConfig.model_validate(data)
+
+
+def run_pretrain(args: argparse.Namespace) -> None:
+    config = load_config(args.config) if args.config else IsMctsConfig()
+    config = _apply_config_overrides(config, args.config_overrides)
+    game_config = config.rules.to_lost_cities_config(seed=config.run.seed)
+    device = torch.device(args.device)
+
+    probe = GameState.new_game(game_config, seed=config.run.seed)
+    in_dim = input_dim(probe, config.encoding)
+    network = AlphaZeroNet.from_config(in_dim, probe.action_size, config).to(device)
+    print(
+        f"network: input_dim={in_dim} action_size={probe.action_size} "
+        f"hidden={config.network.hidden_size} layers={config.network.num_layers}",
+        flush=True,
+    )
+
+    print(f"rolling out {args.games} games of {args.bot} vs {args.bot}...", flush=True)
+    t0 = time.perf_counter()
+    samples = _collect_samples(
+        game_config,
+        config.encoding,
+        n_games=args.games,
+        bot_name=args.bot,
+        seed=args.seed,
+    )
+    rollout_secs = time.perf_counter() - t0
+    print(
+        f"collected {len(samples)} decisions from {args.games} games "
+        f"in {rollout_secs:.1f}s ({len(samples) / args.games:.1f} decisions/game)",
+        flush=True,
+    )
+
+    print(
+        f"training: {args.epochs} epochs, batch={args.batch_size}, lr={args.lr}, "
+        f"value_weight={args.value_loss_weight}",
+        flush=True,
+    )
+    t1 = time.perf_counter()
+    optimizer = _train_supervised(
+        network,
+        samples,
+        device,
+        epochs=args.epochs,
+        batch_size=args.batch_size,
+        lr=args.lr,
+        weight_decay=args.weight_decay,
+        grad_clip=args.grad_clip,
+        value_loss_weight=args.value_loss_weight,
+    )
+    train_secs = time.perf_counter() - t1
+    print(f"training done in {train_secs:.1f}s", flush=True)
+
+    out_path = Path(args.out)
+    out_path.parent.mkdir(parents=True, exist_ok=True)
+    payload = {
+        "config": config.to_dict(),
+        "game_config": game_config.to_snapshot(),
+        "iteration": 0,
+        "network": network.state_dict(),
+        "optimizer": optimizer.state_dict(),
+        "metrics": {
+            "pretrain/games": args.games,
+            "pretrain/samples": len(samples),
+            "pretrain/epochs": args.epochs,
+            "pretrain/bot": args.bot,
+        },
+    }
+    torch.save(payload, out_path)
+    print(f"saved pretrained checkpoint to {out_path}", flush=True)

src/coolrl_lost_cities/games/classic/ismcts/trainer.py

@@ -83,6 +83,51 @@ class IsMctsTrainer:
         self.metrics_path = self.run_dir / "metrics.jsonl"
         self.rng = random.Random(config.run.seed)
 
+        # Optional KL anchor: load a frozen reference network whose policy we
+        # use to regularize updates (KL(current || reference) added to loss).
+        # Prevents drift from a BC pretrained baseline during self-play.
+        self.kl_anchor_ref: AlphaZeroNet | None = None
+        self.kl_anchor_beta = float(config.training.kl_anchor_beta)
+        if config.training.kl_anchor_ckpt and self.kl_anchor_beta > 0.0:
+            ref_path = Path(config.training.kl_anchor_ckpt)
+            if not ref_path.exists():
+                raise FileNotFoundError(f"kl_anchor_ckpt not found: {ref_path}")
+            ref_payload = torch.load(ref_path, map_location=self.device, weights_only=False)
+            self.kl_anchor_ref = AlphaZeroNet.from_config(
+                self.input_dim, self.action_size, config
+            ).to(self.device)
+            self.kl_anchor_ref.load_state_dict(ref_payload["network"])
+            self.kl_anchor_ref.eval()
+            for p in self.kl_anchor_ref.parameters():
+                p.requires_grad = False
+            print(
+                f"[trainer] KL anchor: ref={ref_path} beta={self.kl_anchor_beta}",
+                flush=True,
+            )
+        # Optional mirror-descent reference policy: blend MCTS visit dist
+        # with this reference in log space before computing CE loss.
+        self.md_target_ref: AlphaZeroNet | None = None
+        if config.training.md_target_ref_ckpt:
+            md_ref_path = Path(config.training.md_target_ref_ckpt)
+            if not md_ref_path.exists():
+                raise FileNotFoundError(f"md_target_ref_ckpt not found: {md_ref_path}")
+            md_payload = torch.load(md_ref_path, map_location=self.device, weights_only=False)
+            self.md_target_ref = AlphaZeroNet.from_config(
+                self.input_dim, self.action_size, config
+            ).to(self.device)
+            self.md_target_ref.load_state_dict(md_payload["network"])
+            self.md_target_ref.eval()
+            for p in self.md_target_ref.parameters():
+                p.requires_grad = False
+            print(
+                f"[trainer] mirror-descent ref={md_ref_path} "
+                f"alpha {config.training.md_target_alpha_start} -> "
+                f"{config.training.md_target_alpha_end} over "
+                f"{config.training.md_target_alpha_iters} iters",
+                flush=True,
+            )
+        self._current_md_alpha = float(config.training.md_target_alpha_start)
+
     def _resolve_device(self, device: torch.device | str) -> torch.device:
         token = str(device)
         if token == "auto":
@@ -116,6 +161,15 @@ class IsMctsTrainer:
         return metrics
 
     def run_iteration(self, iteration: int) -> IterationMetrics:
+        # Update mirror-descent alpha schedule (linear from start to end over alpha_iters).
+        if self.md_target_ref is not None:
+            cfg_t = self.config.training
+            n = max(1, int(cfg_t.md_target_alpha_iters))
+            frac = min(1.0, float(iteration) / float(n))
+            self._current_md_alpha = float(
+                cfg_t.md_target_alpha_start
+                + (cfg_t.md_target_alpha_end - cfg_t.md_target_alpha_start) * frac
+            )
         print(
             f"[iter {iteration}] self-play start (workers={self.config.training.num_workers})",
             flush=True,
@@ -258,10 +312,45 @@ class IsMctsTrainer:
         )
         logits, value_pred = self.network(info, legal)
         log_probs = torch.log_softmax(logits, dim=-1)
-        policy_loss = -(pi * log_probs).sum(dim=-1).mean()
+        # Optional mirror-descent target: blend MCTS visit distribution with
+        # the frozen reference (BC) policy in log space, then train CE to that
+        # blended target. This is the standard regularized policy improvement
+        # operator: pi_target = softmax(alpha * log(pi_mcts) + (1-alpha) * log(pi_ref)).
+        if self.md_target_ref is not None:
+            with torch.no_grad():
+                ref_logits, _ref_value = self.md_target_ref(info, legal)
+                ref_log_probs = torch.log_softmax(ref_logits, dim=-1)
+            alpha = float(self._current_md_alpha)
+            # Clamp pi to avoid log(0); MCTS visit dist already has only legal
+            # actions positive, so this affects illegal actions which the mask
+            # in the network forward already zeroed out via -inf logits.
+            log_pi = torch.log(pi.clamp_min(1.0e-12))
+            mixed = alpha * log_pi + (1.0 - alpha) * ref_log_probs
+            mixed = mixed.masked_fill(~legal, torch.finfo(mixed.dtype).min)
+            pi_target = torch.softmax(mixed, dim=-1)
+            policy_loss = -(pi_target * log_probs).sum(dim=-1).mean()
+        else:
+            policy_loss = -(pi * log_probs).sum(dim=-1).mean()
         v_scale = float(self.network.value_scale)
         value_loss = nn.functional.mse_loss(value_pred / v_scale, value_target / v_scale)
-        loss = policy_loss + value_loss
+        value_weight = float(self.config.training.value_loss_weight)
+        loss = policy_loss + value_weight * value_loss
+
+        # KL anchor: KL(current || reference) over legal actions only.
+        # Encourages current policy to stay close to the reference (BC) policy.
+        kl_anchor_loss = 0.0
+        if self.kl_anchor_ref is not None and self.kl_anchor_beta > 0.0:
+            with torch.no_grad():
+                ref_logits, _ref_value = self.kl_anchor_ref(info, legal)
+                ref_log_probs = torch.log_softmax(ref_logits, dim=-1)
+            # KL(current || ref) = sum_a p_cur(a) * (log p_cur(a) - log p_ref(a))
+            cur_probs = log_probs.exp()
+            legal_f = legal.float()
+            kl_per_action = cur_probs * (log_probs - ref_log_probs) * legal_f
+            kl = kl_per_action.sum(dim=-1).mean()
+            loss = loss + self.kl_anchor_beta * kl
+            kl_anchor_loss = float(kl.item())
+
         self.optimizer.zero_grad(set_to_none=True)
         loss.backward()
         if self.config.optimization.grad_clip > 0:
@@ -270,6 +359,9 @@ class IsMctsTrainer:
                 self.config.optimization.grad_clip,
             )
         self.optimizer.step()
+        # Stash auxiliary loss for the IterationMetrics path; we return only
+        # the three primary scalars for backward compatibility.
+        self._last_kl_anchor_loss = kl_anchor_loss
         return float(policy_loss.item()), float(value_loss.item()), float(loss.item())
 
     def _evaluate(self, iteration: int) -> dict[str, float | int]:
@@ -364,9 +456,15 @@ class IsMctsTrainer:
         with torch.inference_mode():
             _logits, value_pred = self.network(info, legal)
             value_error = nn.functional.mse_loss(value_pred, target)
+            value_rmse = float(value_error.item()) ** 0.5
+            target_abs_mean = float(target.abs().mean().item())
+            target_std = float(target.std().item()) if target.numel() > 1 else 0.0
         return {
             "mcts/avg_visit_entropy": float(np.mean(entropies)) if entropies else 0.0,
             "mcts/value_prediction_error": float(value_error.item()),
+            "mcts/value_rmse": value_rmse,
+            "mcts/v_target_abs_mean": target_abs_mean,
+            "mcts/v_target_std": target_std,
             "mcts/policy_mcts_kl": float(np.mean(policy_kls)) if policy_kls else 0.0,
         }
 

tests/games/classic/ismcts/test_ismcts.py

@@ -165,7 +165,13 @@ def test_cython_sequential_matches_python_sequential_visit_counts() -> None:
         )
 
 
-def test_search_visit_counts_match_with_parallel_simulations() -> None:
+def test_search_visit_counts_invariant_with_parallel_simulations() -> None:
+    # Sequential (parallel=1) and batched (parallel=8) MCTS produce different
+    # visit distributions because virtual_loss within a batch spreads simulations
+    # across actions in ways that pure-sequential search does not. The required
+    # invariants are that both legal-action sets and total visit counts match —
+    # this catches the hidden bug where search() ignored parallel_simulations
+    # and always ran batch=1 internally.
     for n_sims in (8, 32, 128):
         state = GameState.new_game(mini_config(), seed=26)
         dim = input_dim(state)
@@ -182,9 +188,17 @@ def test_search_visit_counts_match_with_parallel_simulations() -> None:
             rng=random.Random(28),
         )
 
-        assert batched.search(state, state.current_player) == sequential.search(
+        seq_visits = sequential.search(state, state.current_player)
-            state, state.current_player
+        bat_visits = batched.search(state, state.current_player)
-        )
+        # Same legal action set
+        assert set(seq_visits.keys()) == set(bat_visits.keys())
+        # Both should run a substantial number of sims (early break on terminal
+        # leaf-as-root can leave a few short, but we should be near n_sims).
+        assert sum(seq_visits.values()) >= n_sims - 2
+        assert sum(bat_visits.values()) >= n_sims - 2
+        # Both bounded by n_sims
+        assert sum(seq_visits.values()) <= n_sims
+        assert sum(bat_visits.values()) <= n_sims
 
 
 def test_search_with_virtual_loss_diversity() -> None: