#!/usr/bin/env python3
"""
Torpedo sea-battle simulations.

Reproduces every number in the two Torpedo blog articles:
  * "How to Win at Battleship: 7 Strategies That Actually Work (With the Math)"
  * "Where to Place Your Ships: A Probability Map of the Sea Battle Board"

What it does
  1. Rebuilds the board rules of the Torpedo iOS game: a 10x10 board, the
     classic 5-ship fleet (5, 4, 3, 3, 2; ships may touch) and the "Classic 10"
     fleet (4, 3, 3, 2, 2, 2, 1, 1, 1, 1; ships may not touch, not even at a
     corner, and the water around a sunk ship is revealed).
  2. Re-implements the five computer levels of Torpedo (Beginner, Easy, Normal,
     Hard, Expert) line by line from the Swift source, plus four textbook
     strategies (random, hunt/target, hunt/target with parity, probability
     density), and plays each against fleets drawn uniformly at random from all
     legal layouts.
  3. Tests whether "edge", "centre" or "clustered" fleets survive longer.
  4. Computes where ships end up on the board (occupancy heatmaps) and the best
     first shot.

Every random draw comes from numpy's PCG64 seeded with SEED (printed at start);
work is split into chunks, each with its own child seed derived from SEED and
the chunk's name, so the output is identical whatever the number of workers.

Requirements: Python 3.9+ and numpy.
Usage:  python3 sim.py              # full run, writes results.json, RESULTS.md, *.svg
        python3 sim.py --quick      # small N, for a smoke test (writes *_quick files)
        python3 sim.py --workers 4
"""

import os

# One BLAS thread per process: the work is parallelised across processes.
for _var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "VECLIB_MAXIMUM_THREADS", "MKL_NUM_THREADS"):
    os.environ.setdefault(_var, "1")

import argparse
import json
import math
import multiprocessing as mp
import platform
import sys
import time
import zlib
from collections import Counter

import numpy as np

SEED = 20261007
SIZE = 10
CELLS_N = SIZE * SIZE
COLS = "ABCDEFGHIJ"

# Cell marks, as in CellMark.swift. 0 = not fired at yet.
OPEN, MISS, HIT, WATER = 0, 1, 2, 3

FLEETS = {
    # Fleet.swift: `classic` is played with ships allowed to touch by default
    # (Rules.swift `shipsMayTouch = true`); `classic10` always keeps ships apart
    # (Fleet.requiresShipsApart, Rules.normalized).
    "classic": {"lengths": [5, 4, 3, 3, 2], "touch": True},
    "classic10": {"lengths": [4, 3, 3, 2, 2, 2, 1, 1, 1, 1], "touch": False},
}


def cell_name(i):
    r, c = divmod(int(i), SIZE)
    return f"{COLS[c]}{r + 1}"


def on_board(r, c):
    return 0 <= r < SIZE and 0 <= c < SIZE


# --------------------------------------------------------------------------
# Geometry
# --------------------------------------------------------------------------

ORTH = [[] for _ in range(CELLS_N)]  # orthogonal neighbours of each cell
DIAG_LIST = [[] for _ in range(CELLS_N)]  # diagonal neighbours
for _r in range(SIZE):
    for _c in range(SIZE):
        for _dr, _dc in ((-1, 0), (1, 0), (0, -1), (0, 1)):
            if on_board(_r + _dr, _c + _dc):
                ORTH[_r * SIZE + _c].append((_r + _dr) * SIZE + _c + _dc)
        for _dr, _dc in ((-1, -1), (-1, 1), (1, -1), (1, 1)):
            if on_board(_r + _dr, _c + _dc):
                DIAG_LIST[_r * SIZE + _c].append((_r + _dr) * SIZE + _c + _dc)

DIAG = np.zeros((CELLS_N, CELLS_N))  # DIAG[i, j] = 1 if j is a diagonal neighbour of i
HALO8 = np.zeros((CELLS_N, CELLS_N), bool)  # cell itself plus its 8 neighbours
for _i in range(CELLS_N):
    DIAG[_i, DIAG_LIST[_i]] = 1
    HALO8[_i, _i] = True
    HALO8[_i, ORTH[_i] + DIAG_LIST[_i]] = True

ROW = np.arange(CELLS_N) // SIZE
COL = np.arange(CELLS_N) % SIZE
PARITY2 = (ROW + COL) % 2 == 0
PATTERN = {s: (ROW + COL) % s == 0 for s in range(1, 6)}  # Normal's search pattern


class Placements:
    """Every way one ship of length L fits on the empty board.

    Order follows the Swift loops: Orientation.allCases (horizontal, vertical),
    then Coordinate.all (row-major). With swift=True a 1-square ship is listed
    twice per cell (once per orientation) exactly as ComputerPlayer enumerates
    it; with swift=False duplicates are dropped.
    """

    def __init__(self, length, swift=False):
        cells = []
        for horizontal in (True, False):
            for r in range(SIZE):
                for c in range(SIZE):
                    cs = [(r, c + k) if horizontal else (r + k, c) for k in range(length)]
                    if all(on_board(*x) for x in cs):
                        cells.append([x[0] * SIZE + x[1] for x in cs])
        if length == 1 and not swift:
            cells = cells[: len(cells) // 2]
        self.length = length
        self.cells = np.array(cells, dtype=np.int64)  # (P, L)
        P = len(cells)
        self.count = P
        mb = np.zeros((P, CELLS_N), bool)
        mb[np.arange(P)[:, None], self.cells] = True
        self.mb = mb  # ship cells
        self.m = mb.astype(np.float64)
        self.mt32 = mb.T.astype(np.float32)  # (100, P) for batched matmuls
        ring = np.zeros((P, CELLS_N), bool)  # orthogonal neighbours outside the ship
        halo = np.zeros((P, CELLS_N), bool)  # ship cells plus all 8-neighbours
        for p, cs in enumerate(cells):
            for x in cs:
                ring[p, ORTH[x]] = True
                halo[p] |= HALO8[x]
            ring[p, cs] = False
        self.ringb = ring
        self.ring = ring.astype(np.float64)
        self.halob = halo


GEO = {L: Placements(L) for L in range(1, 6)}
GEO_SWIFT = {L: (Placements(L, swift=True) if L == 1 else GEO[L]) for L in range(1, 6)}


# --------------------------------------------------------------------------
# A battle: hidden fleet + public view
# --------------------------------------------------------------------------

class Game:
    """One side's waters.

    Strategies may read only the public fields, which mirror PublicBoard.swift:
      marks   - per cell OPEN / MISS / HIT / WATER
      unres   - hits on ships still afloat (PublicBoard.unresolvedHits)
      afloat  - lengths of ships still afloat, in fleet order (classes, not places)
      touch   - whether ships may touch
      shots   - shots fired so far (misses + hits; revealed water not counted)
    Sunk ships are announced with their position (Torpedo shows every wreck),
    so their squares stop counting as unresolved hits. Everything prefixed with
    an underscore is hidden information that no strategy touches.
    """

    __slots__ = ("marks", "unres", "afloat", "touch", "shots", "_ids", "_owner", "_cells", "_left")

    def __init__(self, lengths, cells_list, touch):
        self.marks = np.zeros(CELLS_N, np.int8)
        self.unres = np.zeros(CELLS_N, bool)
        self.afloat = list(lengths)
        self.touch = touch
        self.shots = 0
        self._ids = list(range(len(lengths)))
        self._owner = np.full(CELLS_N, -1, np.int16)
        for i, cs in enumerate(cells_list):
            assert (self._owner[cs] == -1).all(), "overlapping ships"
            self._owner[cs] = i
        self._cells = cells_list
        self._left = list(lengths)

    def fire(self, c):
        """Board.fire(at:) in Board.swift."""
        assert self.marks[c] == OPEN, f"strategy fired at a known square {cell_name(c)}"
        self.shots += 1
        o = self._owner[c]
        if o < 0:
            self.marks[c] = MISS
            return 0
        self.marks[c] = HIT
        self.unres[c] = True
        self._left[o] -= 1
        if self._left[o]:
            return 1
        cs = self._cells[o]
        self.unres[cs] = False
        k = self._ids.index(o)
        del self._ids[k]
        del self.afloat[k]
        if not self.touch:  # Board.marksWaterAroundWrecks
            halo = HALO8[cs].any(0)
            self.marks[halo & (self.marks == OPEN)] = WATER
        return 2


def play(strategy, lengths, cells_list, touch, rng, ctx):
    g = Game(lengths, cells_list, touch)
    while g.afloat:
        g.fire(int(strategy(g, rng, ctx)))
        if g.shots > CELLS_N:
            raise RuntimeError("game did not end")
    return g.shots


def pick(rng, seq):
    return seq[int(rng.integers(len(seq)))]


# --------------------------------------------------------------------------
# Reference strategies
# --------------------------------------------------------------------------

def strat_random(g, rng, ctx):
    """Fires at a uniformly random square it hasn't fired at."""
    return pick(rng, np.flatnonzero(g.marks == OPEN))


def _hunt_target(g, rng, parity):
    # Target: probe a random untried neighbour of any hit on a ship still afloat.
    if g.unres.any():
        nb = sorted({n for h in np.flatnonzero(g.unres) for n in ORTH[h] if g.marks[n] == OPEN})
        if nb:
            return pick(rng, nb)
    # Hunt: random square; with parity, only the "black" squares of a
    # checkerboard while every ship afloat is at least 2 long.
    open_ = g.marks == OPEN
    if parity and min(g.afloat) >= 2 and (open_ & PARITY2).any():
        open_ = open_ & PARITY2
    return pick(rng, np.flatnonzero(open_))


def strat_hunt_target(g, rng, ctx):
    return _hunt_target(g, rng, parity=False)


def strat_hunt_target_parity(g, rng, ctx):
    return _hunt_target(g, rng, parity=True)


def strat_density(g, rng, ctx):
    """Probability density (the DataGenetics / Nick Berry approach).

    Every turn, for every ship still afloat, enumerate every placement that is
    consistent with what is known (no miss, no revealed water, no wreck inside
    it; with ships apart, no square diagonal to a live hit and not running
    alongside a live hit it doesn't include). While there are live hits, only
    placements through at least one of them count, weighted by how many they
    cover. Fire at the square covered by the most placements (ties at random).
    """
    hits = g.unres
    open_ = g.marks == OPEN
    anyhit = bool(hits.any())
    hf = hits.astype(np.float64)
    blocked = ~open_ & ~hits
    if not g.touch and anyhit:
        blocked |= (DIAG @ hf) > 0
    bf = blocked.astype(np.float64)
    score = np.zeros(CELLS_N)
    for L, n in Counter(g.afloat).items():
        G = GEO[L]
        valid = (G.m @ bf) == 0
        if anyhit:
            cov = G.m @ hf
            valid &= cov > 0
            if not g.touch:
                valid &= (G.ring @ hf) == 0
            w = np.where(valid, cov, 0.0)
        else:
            w = valid.astype(np.float64)
        score += n * (w @ G.m)
    score[~open_] = 0
    best = score.max()
    if best <= 0:
        return pick(rng, np.flatnonzero(open_))
    return pick(rng, np.flatnonzero(score == best))


# --------------------------------------------------------------------------
# Torpedo's computer player, from ComputerPlayer.swift
# --------------------------------------------------------------------------

def torpedo_random_target(g, rng, follow_up_chance):
    """ComputerPlayer.randomTarget: mostly random, following up on a live hit
    `follow_up_chance` times in ten (Beginner 3, Easy 6)."""
    pend = np.flatnonzero(g.unres)
    if len(pend) and rng.integers(10) < follow_up_chance:
        nb = [n for h in pend for n in ORTH[h] if g.marks[n] == OPEN]  # duplicates kept, as in Swift
        if nb:
            return pick(rng, nb)
    return pick(rng, np.flatnonzero(g.marks == OPEN))


def torpedo_follow_up_candidates(g, pend):
    """ComputerPlayer.followUpCandidates: extend any line of 2+ hits from both
    ends; otherwise every untried orthogonal neighbour of a live hit."""
    if not len(pend):
        return []
    hits = set(pend.tolist())
    cands = []
    for h in pend:
        r, c = divmod(int(h), SIZE)
        for dr, dc in ((0, 1), (1, 0)):
            if not (on_board(r + dr, c + dc) and (r + dr) * SIZE + c + dc in hits):
                continue
            fr, fc = r, c
            while on_board(fr, fc) and fr * SIZE + fc in hits:
                fr, fc = fr + dr, fc + dc
            br, bc = r, c
            while on_board(br, bc) and br * SIZE + bc in hits:
                br, bc = br - dr, bc - dc
            for x, y in ((fr, fc), (br, bc)):
                if on_board(x, y) and g.marks[x * SIZE + y] == OPEN:
                    cands.append(x * SIZE + y)
    if cands:
        return cands
    return [n for h in pend for n in ORTH[h] if g.marks[n] == OPEN]


def torpedo_hunt_and_target(g, rng):
    """ComputerPlayer.huntAndTarget (Normal, and the fallback of Hard/Expert)."""
    pend = np.flatnonzero(g.unres)
    cands = torpedo_follow_up_candidates(g, pend)
    if cands:
        return pick(rng, cands)
    open_ = g.marks == OPEN
    spacing = min(g.afloat)
    # "A gap shorter than the smallest ship afloat can't hide one": a square
    # whose open run across or down is >= spacing is exactly a square covered
    # by some all-open window of that length.
    if spacing == 1:
        fits = open_.copy()
    else:
        G = GEO[spacing]
        valid = (G.m @ (~open_).astype(np.float64)) == 0
        fits = ((valid.astype(np.float64) @ G.m) > 0) & open_
    cand = fits if fits.any() else open_
    pattern = cand & PATTERN[spacing]
    return pick(rng, np.flatnonzero(pattern if pattern.any() else cand))


OPENING_SPREAD_SHOTS = 8  # ComputerPlayer.openingSpreadShots


def torpedo_probability_target(g, rng, habits, spreads_opening):
    """ComputerPlayer.probabilityTarget (Hard and Expert). Returns None when no
    placement fits, so the caller falls back to hunt-and-target."""
    hits = g.unres
    open_ = g.marks == OPEN
    anyhit = bool(hits.any())
    hf = hits.astype(np.float64)
    blocked = ~open_ & ~hits
    if not g.touch:  # knownWater: diagonal neighbours of live hits
        blocked |= (DIAG @ hf) > 0
    bf = blocked.astype(np.float64)
    score = np.zeros(CELLS_N)
    for L, n in Counter(g.afloat).items():
        G = GEO_SWIFT[L]
        valid = (G.m @ bf) == 0
        if anyhit:
            cov = G.m @ hf
            valid &= cov > 0
            if not g.touch:
                valid &= (G.ring @ hf) == 0
            w = np.where(valid, cov * cov, 0.0)  # weight = covered * covered
        else:
            w = valid.astype(np.float64)
        score += n * (w @ G.m)
    score[~open_] = 0
    present = score > 0
    if not present.any():
        return None
    weighted = score
    if not anyhit and habits is not None:
        weighted = score * (1 + 0.6 * habits)
    if spreads_opening and not anyhit and g.shots < OPENING_SPREAD_SHOTS:
        idx = np.flatnonzero(present)
        order = idx[np.lexsort((idx, -weighted[idx]))][:5]  # by weight desc, then (row, column)
        top = weighted[order]
        roll = rng.random() * top.sum()
        for cell, wgt in zip(order, top):
            roll -= wgt
            if roll < 0:
                return int(cell)
    best = weighted[present].max()
    return pick(rng, np.flatnonzero(present & (weighted == best)))


def strat_beginner(g, rng, ctx):
    return torpedo_random_target(g, rng, 3)


def strat_easy(g, rng, ctx):
    return torpedo_random_target(g, rng, 6)


def strat_normal(g, rng, ctx):
    return torpedo_hunt_and_target(g, rng)


def strat_hard(g, rng, ctx):
    t = torpedo_probability_target(g, rng, None, spreads_opening=True)
    return t if t is not None else torpedo_hunt_and_target(g, rng)


def strat_expert(g, rng, ctx):
    t = torpedo_probability_target(g, rng, ctx.get("habits"), spreads_opening=False)
    return t if t is not None else torpedo_hunt_and_target(g, rng)


STRATEGIES = {
    # key: (function, label, group)
    "random": (strat_random, "Random shots", "reference"),
    "hunt_target": (strat_hunt_target, "Hunt/target", "reference"),
    "hunt_target_parity": (strat_hunt_target_parity, "Hunt/target + parity", "reference"),
    "density": (strat_density, "Probability density", "reference"),
    "torpedo_beginner": (strat_beginner, "Torpedo Beginner", "torpedo"),
    "torpedo_easy": (strat_easy, "Torpedo Easy", "torpedo"),
    "torpedo_normal": (strat_normal, "Torpedo Normal", "torpedo"),
    "torpedo_hard": (strat_hard, "Torpedo Hard", "torpedo"),
    "torpedo_expert": (strat_expert, "Torpedo Expert (no history)", "torpedo"),
    "torpedo_expert_learned": (strat_expert, "Torpedo Expert (after 20 games vs. this player)", "torpedo"),
}
HISTORY_GAMES = 20  # games the "learned" Expert has seen the player's fleet in


# --------------------------------------------------------------------------
# Fleet layout samplers
# --------------------------------------------------------------------------

def _cell_mask(L, sel):
    return GEO[L].mb[sel]


def _connected(lengths, idx):
    """Rows where the ships form one orthogonally connected group (each ship
    shares an edge with another, and they all link up)."""
    B, k = idx.shape
    cellm = [_cell_mask(L, idx[:, j]) for j, L in enumerate(lengths)]
    ringm = [GEO[L].ringb[idx[:, j]] for j, L in enumerate(lengths)]
    adj = np.zeros((B, k, k), np.int8)
    for i in range(k):
        for j in range(i + 1, k):
            a = (ringm[i] & cellm[j]).any(1)
            adj[:, i, j] = a
            adj[:, j, i] = a
    reach = np.zeros((B, k), np.int8)
    reach[:, 0] = 1
    for _ in range(k - 1):
        reach = np.maximum(reach, (np.einsum("bi,bij->bj", reach, adj) > 0).astype(np.int8))
    return reach.all(1)


def _style_allowed(style, L):
    """Per-ship placement filters for the placement experiment."""
    G = GEO[L]
    cr, cc = ROW[G.cells], COL[G.cells]
    if style == "uniform" or style == "cluster":
        return np.arange(G.count)
    if style == "edge_ring":  # lies flat along an edge: every square on the outer ring
        ok = ((cr == 0) | (cr == 9) | (cc == 0) | (cc == 9)).all(1)
    elif style == "edge_touch":  # at least one square on the outer ring
        ok = ((cr == 0) | (cr == 9) | (cc == 0) | (cc == 9)).any(1)
    elif style == "centre":  # every square inside the central 6x6 (C3-H8)
        ok = ((cr >= 2) & (cr <= 7) & (cc >= 2) & (cc <= 7)).all(1)
    else:
        raise ValueError(style)
    return np.flatnonzero(ok)


def sample_reject(rng, lengths, n, style="uniform", batch=100_000):
    """Exactly uniform layouts (ships may touch, must not overlap) subject to the
    style's per-ship filter, by rejection: draw each ship independently and
    uniformly among its allowed placements, keep draws with no overlap (and, for
    "cluster", one connected group). Returns placement indices, shape (n, k)."""
    allowed = [_style_allowed(style, L) for L in lengths]
    out, got, proposed = [], 0, 0
    rows = np.arange(batch)[:, None]
    while got < n:
        idx = np.stack([al[rng.integers(len(al), size=batch)] for al in allowed], 1)
        occ = np.zeros((batch, CELLS_N), np.int8)
        for j, L in enumerate(lengths):
            occ[rows, GEO[L].cells[idx[:, j]]] += 1
        ok = occ.max(1) <= 1
        if style == "cluster":
            sub = np.flatnonzero(ok)
            ok[sub] = _connected(lengths, idx[sub])
        proposed += batch
        acc = idx[ok]
        out.append(acc)
        got += len(acc)
    return np.concatenate(out)[:n], got / proposed


def sample_sequential(rng, lengths, touch, n, batch=50_000):
    """Board.random(for:) in Board.swift: place ships in fleet order, each
    uniformly among the spots legal given the ships already placed, starting
    over on a dead end. Also returns log importance weights log(prod n_i),
    which turn these draws into an estimate under the uniform distribution."""
    out, logws, got = [], [], 0
    while got < n:
        H = np.zeros((batch, CELLS_N), np.float32)
        idx = np.zeros((batch, len(lengths)), np.int64)
        logw = np.zeros(batch)
        alive = np.ones(batch, bool)
        for j, L in enumerate(lengths):
            G = GEO[L]
            valid = ((H > 0).astype(np.float32) @ G.mt32) == 0
            cnt = valid.sum(1)
            alive &= cnt > 0
            logw += np.log(np.maximum(cnt, 1))
            u = rng.random(valid.shape, dtype=np.float32)
            u[~valid] = -1
            ch = u.argmax(1)
            idx[:, j] = ch
            H += (G.mb if touch else G.halob)[ch]
        out.append(idx[alive])
        logws.append(logw[alive])
        got += int(alive.sum())
    return np.concatenate(out)[:n], np.concatenate(logws)[:n]


def sample_gibbs(rng, lengths, touch, n, sweeps):
    """Uniform layouts via Gibbs sampling: start from sequential draws, then
    repeatedly move each ship to a uniformly random spot that is legal given the
    others. Uniform over legal layouts is the stationary distribution."""
    idx, _ = sample_sequential(rng, lengths, touch, n)
    B = len(idx)
    H = np.zeros((B, CELLS_N), np.int16)
    blockers = [(GEO[L].mb if touch else GEO[L].halob) for L in lengths]
    for j in range(len(lengths)):
        H += blockers[j][idx[:, j]]
    for _ in range(sweeps):
        for j, L in enumerate(lengths):
            H -= blockers[j][idx[:, j]]
            valid = ((H > 0).astype(np.float32) @ GEO[L].mt32) == 0
            u = rng.random(valid.shape, dtype=np.float32)
            u[~valid] = -1
            ch = u.argmax(1)
            idx[:, j] = ch
            H += blockers[j][ch]
    return idx


def occupancy(lengths, idx, weights=None):
    """Per-cell share of layouts with a ship on the cell."""
    tot = np.zeros(CELLS_N)
    for j, L in enumerate(lengths):
        cells = GEO[L].cells[idx[:, j]]  # (n, L)
        w = None if weights is None else np.repeat(weights, L)
        tot += np.bincount(cells.ravel(), weights=w, minlength=CELLS_N)
    denom = len(idx) if weights is None else weights.sum()
    return tot / denom


def sample_layouts(rng, fleet, style, n):
    f = FLEETS[fleet]
    if style == "sequential":  # Torpedo's own placer, Board.random(for:)
        return sample_sequential(rng, f["lengths"], f["touch"], n)[0]
    if fleet == "classic":
        return sample_reject(rng, f["lengths"], n, style)[0]
    assert style == "uniform"
    return sample_gibbs(rng, f["lengths"], f["touch"], n, GIBBS_SWEEPS)


GIBBS_SWEEPS = 100


# --------------------------------------------------------------------------
# Parallel tasks
# --------------------------------------------------------------------------

def child_rng(name, chunk):
    ss = np.random.SeedSequence(SEED, spawn_key=(zlib.crc32(name.encode()), chunk))
    return np.random.Generator(np.random.PCG64(ss))


def task_games(args):
    """Plays `n` games of one strategy against one kind of fleet."""
    name, fleet, style, strategy, n, chunk = args
    rng = child_rng(name, chunk)
    f = FLEETS[fleet]
    lengths = f["lengths"]
    learned = strategy == "torpedo_expert_learned"
    per_game = 1 + (HISTORY_GAMES if learned else 0)
    layouts = sample_layouts(rng, fleet, style, n * per_game)
    fn = STRATEGIES[strategy][0]
    shots = np.zeros(n, np.int16)
    t0 = time.perf_counter()
    for gi in range(n):
        block = layouts[gi * per_game:(gi + 1) * per_game]
        ctx = {}
        if learned:  # GameModel.habits: share of past games with a ship on each cell
            ctx["habits"] = occupancy(lengths, block[1:]) if HISTORY_GAMES else None
        cells_list = [GEO[L].cells[block[0, j]] for j, L in enumerate(lengths)]
        shots[gi] = play(fn, lengths, cells_list, f["touch"], rng, ctx)
    return name, chunk, shots, time.perf_counter() - t0


def task_heat(args):
    """Draws layouts for the heatmaps; returns per-cell counts."""
    name, kind, n, chunk = args
    rng = child_rng(name, chunk)
    if kind == "classic_uniform":
        idx, acc = sample_reject(rng, FLEETS["classic"]["lengths"], n)
        return name, chunk, {"occ_sum": occupancy(FLEETS["classic"]["lengths"], idx) * n, "n": n, "acceptance": acc}
    if kind in ("classic_sequential", "ten_sequential"):
        fleet = "classic" if kind.startswith("classic") else "classic10"
        f = FLEETS[fleet]
        idx, logw = sample_sequential(rng, f["lengths"], f["touch"], n)
        return name, chunk, {
            "occ_sum": occupancy(f["lengths"], idx) * n, "n": n,
            # importance-weighted sums for the uniform distribution (log-sum-exp safe)
            "logw": logw,
            "occ_w": [np.bincount(GEO[L].cells[idx[:, j]].ravel(), weights=np.repeat(np.exp(logw - 40.0), L), minlength=CELLS_N)
                      for j, L in enumerate(f["lengths"])],
        }
    if kind == "ten_gibbs":
        f = FLEETS["classic10"]
        idx = sample_gibbs(rng, f["lengths"], f["touch"], n, GIBBS_SWEEPS)
        return name, chunk, {"occ_sum": occupancy(f["lengths"], idx) * n, "n": n}
    if kind == "ten_gibbs_short":  # mixing check: far fewer sweeps
        f = FLEETS["classic10"]
        idx = sample_gibbs(rng, f["lengths"], f["touch"], n, 5)
        return name, chunk, {"occ_sum": occupancy(f["lengths"], idx) * n, "n": n}
    if kind == "classic_cluster_rate":
        _, acc = sample_reject(rng, FLEETS["classic"]["lengths"], n, "cluster")
        return name, chunk, {"acceptance": acc}
    raise ValueError(kind)


# --------------------------------------------------------------------------
# Statistics
# --------------------------------------------------------------------------

def summarize(shots):
    a = np.asarray(shots, dtype=np.float64)
    n = len(a)
    hist = np.bincount(a.astype(int), minlength=101)
    return {
        "n": n,
        "mean": round(float(a.mean()), 2),
        "sd": round(float(a.std(ddof=1)), 2),
        "se_mean": round(float(a.std(ddof=1) / math.sqrt(n)), 3),
        "median": float(np.median(a)),
        "p10": float(np.percentile(a, 10)),
        "p90": float(np.percentile(a, 90)),
        "min": int(a.min()),
        "max": int(a.max()),
        "pct_within_50": round(100 * float((a <= 50).mean()), 2),
        "histogram_shots": {str(i): int(c) for i, c in enumerate(hist) if c},
    }


def heat_summary(p, label):
    p = np.asarray(p, dtype=np.float64)
    order = np.argsort(-p, kind="stable")
    centre = [44, 45, 54, 55]  # E5 F5 E6 F6
    corners = [0, 9, 90, 99]  # A1 J1 A10 J10
    hmax, hmin = p.max(), p.min()
    return {
        "label": label,
        "grid_rows_1_to_10_cols_A_to_J": [[round(float(p[r * SIZE + c]), 5) for c in range(SIZE)] for r in range(SIZE)],
        "hottest": [{"cell": cell_name(i), "value": round(float(p[i]), 5)} for i in np.flatnonzero(np.isclose(p, hmax, rtol=0, atol=1e-12))],
        "coldest": [{"cell": cell_name(i), "value": round(float(p[i]), 5)} for i in np.flatnonzero(np.isclose(p, hmin, rtol=0, atol=1e-12))],
        "top5": [{"cell": cell_name(i), "value": round(float(p[i]), 5)} for i in order[:5]],
        "bottom5": [{"cell": cell_name(i), "value": round(float(p[i]), 5)} for i in order[::-1][:5]],
        "centre_mean_E5_F5_E6_F6": round(float(p[centre].mean()), 5),
        "corner_mean_A1_J1_A10_J10": round(float(p[corners].mean()), 5),
        "centre_to_corner_ratio": round(float(p[centre].mean() / p[corners].mean()), 3),
        "max_to_min_ratio": round(float(hmax / hmin), 3),
        "board_mean": round(float(p.mean()), 5),
    }


# --------------------------------------------------------------------------
# SVG output
# --------------------------------------------------------------------------

NAVY, PANEL, TEXT, TEXT2 = "#0B1A2E", "#12253D", "#E8EEF6", "#9FB2C8"
STOPS = ["#12304A", "#2E7FA0", "#4DDEA3"]
FONT = "'Barlow Condensed', 'Arial Narrow', sans-serif"


def _hex(h):
    return [int(h[i:i + 2], 16) for i in (1, 3, 5)]


def scale_colour(t):
    t = min(max(t, 0.0), 1.0)
    a, b, u = (STOPS[0], STOPS[1], t / 0.5) if t <= 0.5 else (STOPS[1], STOPS[2], (t - 0.5) / 0.5)
    ca, cb = _hex(a), _hex(b)
    return "#" + "".join(f"{round(x + (y - x) * u):02X}" for x, y in zip(ca, cb))


def _lum(h):
    r, g, b = (x / 255 for x in _hex(h))
    return 0.2126 * r + 0.7152 * g + 0.0722 * b


def _esc(s):
    return s.replace("&", "&amp;").replace("<", "&lt;").replace(">", "&gt;")


def _contrast(a, b):
    def rl(h):
        cs = [x / 255 for x in _hex(h)]
        cs = [x / 12.92 if x <= 0.03928 else ((x + 0.055) / 1.055) ** 2.4 for x in cs]
        return 0.2126 * cs[0] + 0.7152 * cs[1] + 0.0722 * cs[2]
    la, lb = sorted((rl(a), rl(b)), reverse=True)
    return (la + 0.05) / (lb + 0.05)


def heatmap_svg(p, title, subtitle, desc, note, n_samples):
    p = np.asarray(p)
    lo, hi = float(p.min()), float(p.max())
    cell, gap = 46, 3
    left, top = 46, 112
    grid = SIZE * cell + (SIZE - 1) * gap
    W = left + grid + 28
    H = top + grid + 96
    # Outline the best first shot(s): every square statistically tied with the
    # hottest one (within 2 standard errors of the difference).
    tie = 2 * math.sqrt(2) * math.sqrt(hi * (1 - hi) / n_samples)
    hot = set(np.flatnonzero(p >= hi - tie).tolist())
    o = [
        f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" width="{W}" height="{H}" role="img" aria-labelledby="t d" font-family="{FONT}">',
        f'<title id="t">{_esc(title)}</title>',
        f'<desc id="d">{_esc(desc)}</desc>',
        f'<rect width="{W}" height="{H}" fill="{NAVY}"/>',
        f'<rect x="12" y="12" width="{W - 24}" height="{H - 24}" rx="14" fill="{PANEL}"/>',
        f'<text x="{left - 18}" y="46" fill="{TEXT}" font-size="24" font-weight="600">{_esc(title)}</text>',
        f'<text x="{left - 18}" y="70" fill="{TEXT2}" font-size="15">{_esc(subtitle)}</text>',
    ]
    for c in range(SIZE):
        x = left + c * (cell + gap) + cell / 2
        o.append(f'<text x="{x:.1f}" y="{top - 10}" fill="{TEXT2}" font-size="15" text-anchor="middle">{COLS[c]}</text>')
    for r in range(SIZE):
        y = top + r * (cell + gap) + cell / 2 + 5
        o.append(f'<text x="{left - 10}" y="{y:.1f}" fill="{TEXT2}" font-size="15" text-anchor="end">{r + 1}</text>')
    for i in range(CELLS_N):
        r, c = divmod(i, SIZE)
        x, y = left + c * (cell + gap), top + r * (cell + gap)
        t = 0.0 if hi == lo else (p[i] - lo) / (hi - lo)
        fill = scale_colour(t)
        ink = NAVY if _contrast(fill, NAVY) > _contrast(fill, TEXT) else TEXT
        stroke = f' stroke="{TEXT}" stroke-width="2"' if i in hot else ""
        o.append(f'<rect x="{x}" y="{y}" width="{cell}" height="{cell}" rx="4" fill="{fill}"{stroke}/>')
        o.append(f'<text x="{x + cell / 2:.1f}" y="{y + cell / 2 + 4:.1f}" fill="{ink}" font-size="12" text-anchor="middle">{100 * p[i]:.1f}%</text>')
    # legend
    ly = top + grid + 30
    lw = 240
    o.append('<defs><linearGradient id="g" x1="0" x2="1" y1="0" y2="0">'
             + "".join(f'<stop offset="{k / 2:.1f}" stop-color="{s}"/>' for k, s in enumerate(STOPS))
             + "</linearGradient></defs>")
    o.append(f'<rect x="{left}" y="{ly}" width="{lw}" height="12" rx="3" fill="url(#g)"/>')
    o.append(f'<text x="{left}" y="{ly + 30}" fill="{TEXT2}" font-size="13">{100 * lo:.1f}% (coldest)</text>')
    o.append(f'<text x="{left + lw}" y="{ly + 30}" fill="{TEXT2}" font-size="13" text-anchor="end">{100 * hi:.1f}% (hottest)</text>')
    o.append(f'<text x="{left + grid}" y="{ly + 10}" fill="{TEXT2}" font-size="13" text-anchor="end">{_esc(note)}</text>')
    o.append(f'<text x="{left + grid}" y="{ly + 30}" fill="{TEXT2}" font-size="13" text-anchor="end">Outlined: best first shot{'s (tied)' if len(hot) > 1 else ''}</text>')
    o.append("</svg>")
    return "\n".join(o) + "\n"


def bars_svg(rows, title, subtitle, desc, note):
    """rows: list of (label, mean, group) sorted best first."""
    label_w, left_pad, bar_max = 250, 28, 420
    row_h, top = 40, 104
    W = left_pad + label_w + bar_max + 110
    H = top + row_h * len(rows) + 86
    xmax = 100.0
    x0 = left_pad + label_w
    o = [
        f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" width="{W}" height="{H}" role="img" aria-labelledby="t d" font-family="{FONT}">',
        f'<title id="t">{_esc(title)}</title>',
        f'<desc id="d">{_esc(desc)}</desc>',
        f'<rect width="{W}" height="{H}" fill="{NAVY}"/>',
        f'<rect x="12" y="12" width="{W - 24}" height="{H - 24}" rx="14" fill="{PANEL}"/>',
        f'<text x="{left_pad}" y="46" fill="{TEXT}" font-size="24" font-weight="600">{_esc(title)}</text>',
        f'<text x="{left_pad}" y="70" fill="{TEXT2}" font-size="15">{_esc(subtitle)}</text>',
    ]
    for v in (0, 25, 50, 75, 100):
        x = x0 + bar_max * v / xmax
        o.append(f'<line x1="{x:.1f}" x2="{x:.1f}" y1="{top - 6}" y2="{top + row_h * len(rows)}" stroke="{NAVY}" stroke-width="1.5"/>')
        o.append(f'<text x="{x:.1f}" y="{top + row_h * len(rows) + 20}" fill="{TEXT2}" font-size="13" text-anchor="middle">{v}</text>')
    o.append(f'<text x="{x0 + bar_max / 2:.1f}" y="{top + row_h * len(rows) + 40}" fill="{TEXT2}" font-size="13" text-anchor="middle">average shots to sink the whole fleet (fewer is better)</text>')
    for k, (label, mean, group) in enumerate(rows):
        y = top + k * row_h
        colour = STOPS[2] if group == "torpedo" else STOPS[1]
        w = bar_max * mean / xmax
        o.append(f'<text x="{x0 - 12}" y="{y + 22}" fill="{TEXT}" font-size="16" text-anchor="end">{_esc(label)}</text>')
        o.append(f'<rect x="{x0}" y="{y + 6}" width="{w:.1f}" height="{row_h - 14}" rx="4" fill="{colour}"/>')
        o.append(f'<text x="{x0 + w + 8:.1f}" y="{y + 22}" fill="{TEXT}" font-size="16" font-weight="600">{mean:.1f}</text>')
    ly = H - 30
    o.append(f'<rect x="{left_pad}" y="{ly - 11}" width="12" height="12" rx="2" fill="{STOPS[1]}"/>')
    o.append(f'<text x="{left_pad + 18}" y="{ly}" fill="{TEXT2}" font-size="13">Textbook strategy</text>')
    o.append(f'<rect x="{left_pad + 140}" y="{ly - 11}" width="12" height="12" rx="2" fill="{STOPS[2]}"/>')
    o.append(f'<text x="{left_pad + 158}" y="{ly}" fill="{TEXT2}" font-size="13">Torpedo computer level</text>')
    o.append(f'<text x="{W - 28}" y="{ly}" fill="{TEXT2}" font-size="13" text-anchor="end">{_esc(note)}</text>')
    o.append("</svg>")
    return "\n".join(o) + "\n"


# --------------------------------------------------------------------------
# Main
# --------------------------------------------------------------------------

MAIN_STRATEGIES = ["random", "hunt_target", "hunt_target_parity", "density",
                   "torpedo_beginner", "torpedo_easy", "torpedo_normal", "torpedo_hard", "torpedo_expert"]
TORPEDO_LEVELS = ["torpedo_beginner", "torpedo_easy", "torpedo_normal", "torpedo_hard", "torpedo_expert"]
PLACEMENT_STYLES = ["uniform", "edge_ring", "edge_touch", "centre", "cluster"]
PLACEMENT_STRATEGIES = ["hunt_target_parity", "density", "torpedo_hard", "torpedo_expert_learned"]
STYLE_LABELS = {
    "uniform": "Uniform random legal layout",
    "edge_ring": "Edge-hugging: every ship lies flat along an edge (all squares on the outer ring)",
    "edge_touch": "Edge-touching: every ship has at least one square on the outer ring",
    "centre": "Centre cluster: every square inside the central 6x6 (C3-H8)",
    "cluster": "All ships touching: one orthogonally connected group of ships",
}


def plan(quick):
    s = 0.02 if quick else 1.0
    n = lambda x: max(40, int(x * s))
    chunk = 250 if not quick else 40
    return {
        "chunk": chunk,
        "main_classic": {k: n(100_000) for k in MAIN_STRATEGIES},
        "main_ten": {k: n(50_000) for k in MAIN_STRATEGIES},
        "placement": n(20_000),
        "validation": n(20_000),
        "heat": {
            "classic_uniform": (n(4_000_000), 200_000 if not quick else 40_000),
            "classic_sequential": (n(2_000_000), 100_000 if not quick else 20_000),
            "ten_sequential": (n(4_000_000), 100_000 if not quick else 20_000),
            "ten_gibbs": (n(1_000_000), 50_000 if not quick else 10_000),
            "ten_gibbs_short": (n(200_000), 50_000 if not quick else 10_000),
            "classic_cluster_rate": (n(20_000), 2_000 if not quick else 400),
        },
    }


def chunks(total, size):
    out, k = [], 0
    while total > 0:
        out.append((k, min(size, total)))
        total -= size
        k += 1
    return out


def empty_board_density(lengths, swift=False):
    geo = GEO_SWIFT if swift else GEO
    tot = np.zeros(CELLS_N)
    for L in lengths:
        tot += geo[L].m.sum(0)
    return tot


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--quick", action="store_true", help="tiny N smoke test; writes *_quick outputs")
    ap.add_argument("--workers", type=int, default=max(1, (os.cpu_count() or 2) - 2))
    ap.add_argument("--out", default=os.path.dirname(os.path.abspath(__file__)))
    args = ap.parse_args()
    P = plan(args.quick)
    suffix = "_quick" if args.quick else ""
    print(f"seed = {SEED}  workers = {args.workers}  quick = {args.quick}", flush=True)
    t_start = time.time()
    timings = {}

    # ---- build task list -------------------------------------------------
    game_tasks = []
    for strat, n in P["main_classic"].items():
        for k, m in chunks(n, P["chunk"]):
            game_tasks.append((f"classic/uniform/{strat}", "classic", "uniform", strat, m, k))
    for strat, n in P["main_ten"].items():
        for k, m in chunks(n, P["chunk"]):
            game_tasks.append((f"classic10/uniform/{strat}", "classic10", "uniform", strat, m, k))
    for style in PLACEMENT_STYLES:
        for strat in PLACEMENT_STRATEGIES:
            for k, m in chunks(P["placement"], P["chunk"]):
                game_tasks.append((f"placement/{style}/{strat}", "classic", style, strat, m, k))
    for fleet in FLEETS:
        for strat in TORPEDO_LEVELS:
            for k, m in chunks(P["validation"], P["chunk"]):
                game_tasks.append((f"validation/{fleet}/{strat}", fleet, "sequential", strat, m, k))
    heat_tasks = []
    for kind, (n, size) in P["heat"].items():
        for k, m in chunks(n, size):
            heat_tasks.append((f"heat/{kind}", kind, m, k))

    # Slow tasks first for better load balance.
    slow = {"density", "torpedo_hard", "torpedo_expert", "torpedo_expert_learned"}
    game_tasks.sort(key=lambda t: (t[3] not in slow, t[1] != "classic"))

    ctx = mp.get_context("spawn")
    game_res, heat_res, cpu = {}, {}, {}
    with ctx.Pool(args.workers) as pool:
        t0 = time.time()
        for name, chunk, payload in pool.imap_unordered(task_heat, heat_tasks):
            heat_res.setdefault(name, {})[chunk] = payload
        timings["heatmaps_wall_s"] = round(time.time() - t0, 1)
        print(f"heatmaps done in {timings['heatmaps_wall_s']} s", flush=True)
        t0 = time.time()
        done = 0
        for name, chunk, shots, secs in pool.imap_unordered(task_games, game_tasks, chunksize=1):
            game_res.setdefault(name, {})[chunk] = shots
            cpu[name] = cpu.get(name, 0.0) + secs
            done += 1
            if done % 50 == 0 or done == len(game_tasks):
                print(f"  games: {done}/{len(game_tasks)} chunks  ({time.time() - t0:.0f} s)", flush=True)
        timings["games_wall_s"] = round(time.time() - t0, 1)

    def shots_of(name):
        parts = game_res[name]
        return np.concatenate([parts[k] for k in sorted(parts)])

    results = {
        "meta": {
            "seed": SEED,
            "rng": "numpy PCG64; child streams SeedSequence(seed, spawn_key=(crc32(task name), chunk index))",
            "quick": args.quick,
            "python": platform.python_version(),
            "numpy": np.__version__,
            "workers": args.workers,
            "chunk_games": P["chunk"],
            "rules": {
                "board": "10x10",
                "classic": "5 ships of 5, 4, 3, 3, 2 squares; ships may touch (Torpedo default) but not overlap",
                "classic10": "10 ships of 4, 3, 3, 2, 2, 2, 1, 1, 1, 1 squares; ships may not touch, even diagonally; water around a sunk ship is revealed",
                "shots": "one shot at a time; every shot counts (Torpedo's default Hit Streak rule gives extra turns but never changes how many shots a side needs)",
                "information": "hit/miss reported; a sunk ship is announced with its class and squares (Torpedo shows every wreck)",
            },
            "gibbs_sweeps": GIBBS_SWEEPS,
            "history_games_for_learned_expert": HISTORY_GAMES,
        },
    }

    # ---- games -------------------------------------------------------------
    for fleet, key in (("classic", "main_classic"), ("classic10", "main_ten")):
        block = {}
        for strat in MAIN_STRATEGIES:
            name = f"{fleet}/uniform/{strat}"
            s = summarize(shots_of(name))
            s["label"] = STRATEGIES[strat][1]
            s["group"] = STRATEGIES[strat][2]
            s["cpu_seconds"] = round(cpu[name], 1)
            s["ms_per_game"] = round(1000 * cpu[name] / s["n"], 2)
            block[strat] = s
        results[f"games_{fleet}"] = block

    plc = {"styles": STYLE_LABELS, "n_per_cell": P["placement"], "table": {}}
    for style in PLACEMENT_STYLES:
        row = {}
        for strat in PLACEMENT_STRATEGIES:
            s = summarize(shots_of(f"placement/{style}/{strat}"))
            s.pop("histogram_shots")
            s["label"] = STRATEGIES[strat][1]
            row[strat] = s
        plc["table"][style] = row
    results["placement_experiment"] = plc

    # Cross-check against Torpedo's real Swift code (swift-check/), both on
    # fleets from Board.random, the app's own placer.
    swift = {}
    sf = os.path.join(os.path.dirname(os.path.abspath(__file__)), "swift-check", "swift_harness.json")
    if os.path.exists(sf):
        with open(sf) as fh:
            swift = json.load(fh)
    val = {"note": "Python port vs. the app's own Swift ComputerPlayer, both against fleets from Board.random (not uniform)", "rows": {}}
    for fleet in FLEETS:
        for strat in TORPEDO_LEVELS:
            s = summarize(shots_of(f"validation/{fleet}/{strat}"))
            row = {"python_n": s["n"], "python_mean": s["mean"], "python_se": s["se_mean"]}
            key = f"{fleet}/{strat.replace('torpedo_', '')}"
            if key in swift:
                sw = swift[key]
                row.update({"swift_n": sw["n"], "swift_mean": round(sw["mean"], 2), "swift_se": round(sw["se"], 3)})
                row["z"] = round((s["mean"] - sw["mean"]) / math.hypot(s["se_mean"], sw["se"]), 2)
            val["rows"][key] = row
    results["validation_vs_swift"] = val

    # ---- heatmaps ----------------------------------------------------------
    def merged(kind):
        parts = heat_res[f"heat/{kind}"]
        return [parts[k] for k in sorted(parts)]

    heat = {}
    cl = merged("classic_uniform")
    n_cl = sum(x["n"] for x in cl)
    p_cl = sum(x["occ_sum"] for x in cl) / n_cl
    heat["classic_uniform"] = heat_summary(p_cl, "Classic 5-ship fleet, uniform random legal layouts (rejection sampling, exact uniform)")
    heat["classic_uniform"]["n_layouts"] = n_cl
    heat["classic_uniform"]["max_standard_error"] = round(float(np.sqrt(p_cl * (1 - p_cl) / n_cl).max()), 6)
    heat["classic_uniform"]["rejection_acceptance_rate"] = round(float(np.mean([x["acceptance"] for x in cl])), 4)

    for kind, fleet in (("classic_sequential", "classic"), ("ten_sequential", "classic10")):
        parts = merged(kind)
        n = sum(x["n"] for x in parts)
        p_seq = sum(x["occ_sum"] for x in parts) / n
        lengths = FLEETS[fleet]["lengths"]
        # importance-weighted estimate of the uniform distribution
        w_all = np.concatenate([np.exp(x["logw"] - 40.0) for x in parts])
        ess = float(w_all.sum() ** 2 / (w_all ** 2).sum())
        occ_w = sum(sum(x["occ_w"]) for x in parts) / w_all.sum()
        heat[kind] = heat_summary(p_seq, f"{fleet}: Torpedo's own random placer (Board.random, ships placed one by one)")
        heat[kind]["n_layouts"] = n
        heat[f"{kind}_reweighted_to_uniform"] = heat_summary(occ_w, f"{fleet}: sequential draws re-weighted to uniform (importance sampling)")
        heat[f"{kind}_reweighted_to_uniform"]["n_layouts"] = n
        heat[f"{kind}_reweighted_to_uniform"]["effective_sample_size"] = round(ess)

    tg = merged("ten_gibbs")
    n_tg = sum(x["n"] for x in tg)
    p_tg = sum(x["occ_sum"] for x in tg) / n_tg
    heat["ten_uniform"] = heat_summary(p_tg, f"Classic 10 fleet (no touching), uniform random legal layouts (Gibbs sampler, {GIBBS_SWEEPS} sweeps)")
    heat["ten_uniform"]["n_layouts"] = n_tg
    heat["ten_uniform"]["max_standard_error_if_independent"] = round(float(np.sqrt(p_tg * (1 - p_tg) / n_tg).max()), 6)
    ts = merged("ten_gibbs_short")
    p_ts = sum(x["occ_sum"] for x in ts) / sum(x["n"] for x in ts)
    p_tw = np.array([v for row in heat["ten_sequential_reweighted_to_uniform"]["grid_rows_1_to_10_cols_A_to_J"] for v in row])
    heat["ten_uniform"]["check_max_abs_diff_vs_5_sweeps"] = round(float(np.abs(p_ts - p_tg).max()), 5)
    heat["ten_uniform"]["check_max_abs_diff_vs_importance_sampling"] = round(float(np.abs(p_tw - p_tg).max()), 5)
    p_cs = np.array([v for row in heat["classic_sequential_reweighted_to_uniform"]["grid_rows_1_to_10_cols_A_to_J"] for v in row])
    heat["classic_uniform"]["check_max_abs_diff_vs_importance_sampling"] = round(float(np.abs(p_cs - p_cl).max()), 5)

    cr = merged("classic_cluster_rate")
    heat["cluster_layout_share"] = {
        "share_of_uniform_layouts_that_form_one_connected_group": round(float(np.mean([x["acceptance"] for x in cr])) / heat["classic_uniform"]["rejection_acceptance_rate"], 5),
        "note": "acceptance of the cluster rejection sampler divided by the plain no-overlap acceptance",
    }

    dens_cl = empty_board_density(FLEETS["classic"]["lengths"])
    dens_ten = empty_board_density(FLEETS["classic10"]["lengths"])
    heat["empty_board_density_classic"] = heat_summary(dens_cl, "Classic fleet: number of single-ship placements covering each square on an empty board (exact count, summed over the 5 ships)")
    heat["empty_board_density_classic"]["total_placements"] = int(sum(GEO[L].count for L in FLEETS["classic"]["lengths"]))
    heat["empty_board_density_ten"] = heat_summary(dens_ten, "Classic 10 fleet: single-ship placements covering each square (exact count, summed over the 10 ships, 1-square ships counted once per square)")
    heat["empty_board_density_carrier"] = heat_summary(GEO[5].m.sum(0), "Carrier (5) alone: placements covering each square")
    heat["empty_board_density_destroyer"] = heat_summary(GEO[2].m.sum(0), "Destroyer (2) alone: placements covering each square")
    results["heatmaps"] = heat

    # best first shot
    best_cl = float(p_cl.max())
    results["best_first_shot"] = {
        "classic": {
            "cells": [cell_name(i) for i in np.flatnonzero(np.isclose(p_cl, best_cl, atol=1e-12))],
            "statistically_tied_within_2_se": [cell_name(i) for i in np.flatnonzero(p_cl >= best_cl - 2 * np.sqrt(best_cl * (1 - best_cl) / n_cl) * math.sqrt(2))],
            "hit_probability": round(best_cl, 4),
            "board_average_hit_probability": round(17 / 100, 4),
            "corner_hit_probability_mean": heat["classic_uniform"]["corner_mean_A1_J1_A10_J10"],
            "density_argmax_cells": heat["empty_board_density_classic"]["hottest"],
        },
        "classic10": {
            "cells": [cell_name(i) for i in np.flatnonzero(np.isclose(p_tg, p_tg.max(), atol=1e-12))],
            "statistically_tied_within_2_se": [cell_name(i) for i in np.flatnonzero(p_tg >= p_tg.max() - 2 * np.sqrt(p_tg.max() * (1 - p_tg.max()) / n_tg) * math.sqrt(2))],
            "hit_probability": round(float(p_tg.max()), 4),
            "board_average_hit_probability": round(20 / 100, 4),
            "corner_hit_probability_mean": heat["ten_uniform"]["corner_mean_A1_J1_A10_J10"],
        },
    }

    timings["total_wall_s"] = round(time.time() - t_start, 1)
    results["meta"]["timings"] = timings

    out = args.out
    with open(os.path.join(out, f"results{suffix}.json"), "w") as fh:
        json.dump(results, fh, indent=1)

    # ---- SVGs ----------------------------------------------------------------
    hc = heat["classic_uniform"]
    with open(os.path.join(out, f"heatmap-classic{suffix}.svg"), "w") as fh:
        fh.write(heatmap_svg(
            p_cl,
            "Where ships end up: classic 5-ship fleet",
            f"Chance each square holds a ship · {n_cl:,} random legal layouts · seed {SEED}",
            (f"10 by 10 heatmap, columns A to J, rows 1 to 10. Probability that each square is occupied when the "
             f"classic fleet (5, 4, 3, 3 and 2 squares, ships may touch) is placed uniformly at random among all legal layouts. "
             f"Hottest: {', '.join(x['cell'] for x in hc['top5'][:2])} at {100 * hc['top5'][0]['value']:.1f}%. "
             f"Coldest: corners at {100 * hc['corner_mean_A1_J1_A10_J10']:.1f}%. Board average 17%."),
            "17 ship squares of 100: average 17%", n_cl))
    ht = heat["ten_uniform"]
    with open(os.path.join(out, f"heatmap-ten-ship{suffix}.svg"), "w") as fh:
        fh.write(heatmap_svg(
            p_tg,
            "Where ships end up: 10-ship fleet, no touching",
            f"Chance each square holds a ship · {n_tg:,} random legal layouts · seed {SEED}",
            (f"10 by 10 heatmap, columns A to J, rows 1 to 10. Probability that each square is occupied when the "
             f"10-ship fleet (4, 3, 3, 2, 2, 2, 1, 1, 1, 1 squares, ships may not touch) is placed uniformly at random among all legal layouts. "
             f"Hottest: edge squares near, not at, the corners ({', '.join(x['cell'] for x in ht['top5'])}) at about {100 * ht['top5'][0]['value']:.1f}%, against {100 * ht['centre_mean_E5_F5_E6_F6']:.1f}% in the centre. "
             f"Coldest: {', '.join(x['cell'] for x in ht['bottom5'][:3])} at about {100 * ht['bottom5'][0]['value']:.1f}%. Board average 20%."),
            "20 ship squares of 100: average 20%", n_tg))
    gc = results["games_classic"]
    rows = sorted(((gc[k]["label"], gc[k]["mean"], gc[k]["group"]) for k in MAIN_STRATEGIES), key=lambda r: r[1])
    with open(os.path.join(out, f"avg-shots-by-strategy{suffix}.svg"), "w") as fh:
        fh.write(bars_svg(
            rows,
            "How many shots to sink the classic fleet?",
            f"Average over {min(gc[k]['n'] for k in MAIN_STRATEGIES):,}+ simulated games per strategy · seed {SEED}",
            "Horizontal bar chart of average shots needed to sink the classic 5-ship fleet: "
            + "; ".join(f"{label} {mean:.1f}" for label, mean, _ in rows) + ".",
            "Classic fleet 5-4-3-3-2, random legal layouts"))

    write_markdown(results, os.path.join(out, f"RESULTS{suffix}.md"))
    print(f"done in {timings['total_wall_s']} s", flush=True)


# --------------------------------------------------------------------------
# RESULTS.md
# --------------------------------------------------------------------------

def _pct(x, d=1):
    return f"{100 * x:.{d}f}%"


def write_markdown(R, path):
    m = R["meta"]
    L = []
    w = L.append
    w("# Torpedo sea-battle simulations: results\n")
    w(f"Seed **{m['seed']}** · numpy {m['numpy']} · Python {m['python']} · {m['workers']} worker processes · "
      f"total run time {m['timings']['total_wall_s']} s (heatmaps {m['timings']['heatmaps_wall_s']} s, games {m['timings']['games_wall_s']} s).")
    w("Every number below is printed by `sim.py` from `results.json`; rerunning `python3 sim.py` reproduces them exactly.\n")

    w("## Rules simulated\n")
    for k, v in m["rules"].items():
        w(f"- **{k}**: {v}")
    w("")

    for fleet, title in (("classic", "Classic 5-ship fleet (5-4-3-3-2, ships may touch)"),
                         ("classic10", "Classic 10 fleet (4-3-3-2-2-2-1-1-1-1, no touching)")):
        g = R[f"games_{fleet}"]
        w(f"## Shots to sink the whole fleet: {title}\n")
        w("Fleets drawn uniformly at random from all legal layouts. Fewer shots is better.\n")
        w("| Strategy | N | Mean (± s.e.) | Median | p10 | p90 | Min | Max | Done in ≤ 50 shots |")
        w("|---|---:|---:|---:|---:|---:|---:|---:|---:|")
        for k in MAIN_STRATEGIES:
            s = g[k]
            w(f"| {s['label']} | {s['n']:,} | {s['mean']:.2f} ± {s['se_mean']:.2f} | {s['median']:g} | {s['p10']:g} | {s['p90']:g} | {s['min']} | {s['max']} | {s['pct_within_50']:.1f}% |")
        w("")

    p = R["placement_experiment"]
    w("## Placement experiment: are edges safer? (classic fleet)\n")
    w(f"N = {p['n_per_cell']:,} games per cell. Each fleet is drawn uniformly at random among the legal layouts that satisfy the style. "
      "The \"learned\" Expert first sees 20 earlier fleets of the same style, as Torpedo's Expert would after 20 games against that player.\n")
    for k, v in p["styles"].items():
        w(f"- **{k}**: {v}")
    w("")
    w("| Fleet style | " + " | ".join(STRATEGIES[s][1] for s in PLACEMENT_STRATEGIES) + " |")
    w("|---|" + "---:|" * len(PLACEMENT_STRATEGIES))
    for style, row in p["table"].items():
        w(f"| {style} | " + " | ".join(f"{row[s]['mean']:.2f} ± {row[s]['se_mean']:.2f}" for s in PLACEMENT_STRATEGIES) + " |")
    w("")
    w("Share of games finished within 50 shots:\n")
    w("| Fleet style | " + " | ".join(STRATEGIES[s][1] for s in PLACEMENT_STRATEGIES) + " |")
    w("|---|" + "---:|" * len(PLACEMENT_STRATEGIES))
    for style, row in p["table"].items():
        w(f"| {style} | " + " | ".join(f"{row[s]['pct_within_50']:.1f}%" for s in PLACEMENT_STRATEGIES) + " |")
    w("")

    v = R["validation_vs_swift"]
    w("## Cross-check against Torpedo's Swift code\n")
    w("`swift-check/` compiles Torpedo's own `ComputerPlayer.target` and board code (unmodified copies) into a command-line "
      "harness. Both it and `sim.py` play each level against fleets from `Board.random` (the app's own placer). "
      "z = difference of means divided by its standard error; |z| < 2 means the port is statistically indistinguishable.\n")
    w("| Fleet / level | Python N | Python mean | Swift N | Swift mean | z |")
    w("|---|---:|---:|---:|---:|---:|")
    for k, r in v["rows"].items():
        w(f"| {k} | {r['python_n']:,} | {r['python_mean']:.2f} ± {r['python_se']:.2f} | "
          + (f"{r['swift_n']:,} | {r['swift_mean']:.2f} ± {r['swift_se']:.2f} | {r['z']:+.2f} |" if "swift_n" in r else "– | – | – |"))
    w("")

    H = R["heatmaps"]
    w("## Heatmaps: where ships end up\n")
    for key in ("classic_uniform", "ten_uniform", "classic_sequential", "ten_sequential",
                "empty_board_density_classic", "empty_board_density_ten"):
        h = H[key]
        is_count = key.startswith("empty")
        fmt = (lambda v: f"{v:g}") if is_count else (lambda v: _pct(v, 2))
        w(f"### {h['label']}\n")
        if "n_layouts" in h:
            w(f"- N = {h['n_layouts']:,} layouts")
        w(f"- Hottest: {', '.join(x['cell'] + ' ' + fmt(x['value']) for x in h['hottest'])}; top 5: {', '.join(x['cell'] + ' ' + fmt(x['value']) for x in h['top5'])}")
        w(f"- Coldest: {', '.join(x['cell'] + ' ' + fmt(x['value']) for x in h['coldest'])}")
        w(f"- Centre (E5, F5, E6, F6) mean {fmt(h['centre_mean_E5_F5_E6_F6'])} vs corner (A1, J1, A10, J10) mean {fmt(h['corner_mean_A1_J1_A10_J10'])}: ratio **{h['centre_to_corner_ratio']}**; hottest/coldest ratio {h['max_to_min_ratio']}")
        for extra in ("max_standard_error", "rejection_acceptance_rate", "check_max_abs_diff_vs_importance_sampling",
                      "max_standard_error_if_independent", "check_max_abs_diff_vs_5_sweeps", "total_placements"):
            if extra in h:
                w(f"- {extra.replace('_', ' ')}: {h[extra]}")
        w("")
        w("| | " + " | ".join(COLS) + " |")
        w("|---|" + "---:|" * SIZE)
        for r, row in enumerate(h["grid_rows_1_to_10_cols_A_to_J"]):
            w(f"| **{r + 1}** | " + " | ".join((f"{v:g}" if is_count else f"{100 * v:.1f}") for v in row) + " |")
        w("")
        if not is_count:
            w("(Values in %.)\n")
    for key in ("classic_sequential_reweighted_to_uniform", "ten_sequential_reweighted_to_uniform"):
        h = H[key]
        w(f"- Cross-check, {h['label']}: hottest {', '.join(x['cell'] + ' ' + _pct(x['value'], 2) for x in h['hottest'])}, effective sample size {h['effective_sample_size']:,} of {h['n_layouts']:,}.")
    c = H["cluster_layout_share"]
    w(f"- Share of uniform random classic layouts in which all five ships form one connected group: {_pct(c['share_of_uniform_layouts_that_form_one_connected_group'], 2)}.\n")

    b = R["best_first_shot"]
    w("## Best first shot\n")
    w(f"- Classic fleet: **{', '.join(b['classic']['cells'])}** with hit probability **{_pct(b['classic']['hit_probability'], 2)}** "
      f"(board average {_pct(b['classic']['board_average_hit_probability'], 0)}, corners {_pct(b['classic']['corner_hit_probability_mean'], 2)}). "
      f"Squares statistically tied with it (within 2 standard errors of the difference): {', '.join(b['classic']['statistically_tied_within_2_se'])}. "
      f"The empty-board placement count (the density hunter's first pick) peaks at {', '.join(x['cell'] for x in b['classic']['density_argmax_cells'])}.")
    w(f"- Classic 10 fleet: **{', '.join(b['classic10']['cells'])}** with hit probability **{_pct(b['classic10']['hit_probability'], 2)}** "
      f"(board average {_pct(b['classic10']['board_average_hit_probability'], 0)}, corners {_pct(b['classic10']['corner_hit_probability_mean'], 2)}). "
      f"Statistically tied: {', '.join(b['classic10']['statistically_tied_within_2_se'])}.\n")

    w(METHODOLOGY)
    w(AI_MODELLING)
    w(OBSERVATIONS)
    with open(path, "w") as fh:
        fh.write("\n".join(L) + "\n")


METHODOLOGY = """## Methodology

The board, fleets and shot rules are rebuilt from Torpedo's Swift model (`app/App/Model`). A ship is a straight
line of squares, horizontal or vertical, fully on the 10x10 board; ships never overlap. In the classic fleet ships
may touch (Torpedo's default); in Classic 10 they may not touch even at a corner, and when a ship sinks the water
around it is marked, exactly as `Board.fire(at:)` does. Every strategy, including Torpedo's own computer levels, sees
only the public view (`PublicBoard.swift`): hits, misses, revealed water, and each sunk ship's class and squares
(Torpedo shows every wreck). Turn order does not matter for these figures: we count the shots one side needs to
sink the other side's whole fleet, which is the same under One Shot, Hit Streak or any other per-shot rule.

Fleets are drawn **uniformly at random among all legal layouts**. For the classic fleet this is exact: each ship
is drawn independently and uniformly among its positions, and draws with an overlap are rejected (also used for the
placement styles, with a per-ship filter, and with a connectivity check for "all ships touching"). For the
no-touch 10-ship fleet rejection is too slow, so we use a Gibbs sampler (start from a sequential draw, then
100 sweeps in which each ship in turn moves to a uniformly random spot that is legal given the others); its heatmap
is checked against a 5-sweep run and against importance-weighted sequential draws (differences listed above).
Torpedo's own auto-placer, `Board.random(for:)`, places ships one at a time and is *not* uniform; its heatmaps are
reported separately. Statistics: mean with standard error, median and 10th/90th percentiles (numpy linear
interpolation), min, max, and the share of games won in at most 50 shots. Heatmap probabilities are per-square
occupancy frequencies over the sampled layouts; the empty-board density maps are exact counts.

The "probability density" reference strategy is the per-ship placement count popularised by Nick Berry
(DataGenetics, 2011): every turn it counts, for each ship still afloat, every placement consistent with the known
misses, revealed water and wrecks, restricted to placements through a live hit while one exists (weighted by the
number of live hits covered), and fires at the square covered by the most placements. It counts placements of each
ship separately, not whole joint fleet layouts.
"""

AI_MODELLING = """## How each Torpedo computer level was modelled

All references are to `app/App/Model/` unless noted. Torpedo's computer aims through `PublicBoard`
(PublicBoard.swift:3-34), so it never sees afloat ships; the fairness tests in
`app/Tests/ComputerFairnessTests.swift:52-98` check that moving a hidden, untouched ship never changes any level's move.
The simulator gives the strategies the same public fields only.

- **Beginner**: `randomTarget(followUpChance: 3)` (ComputerPlayer.swift:31-32, 271-286). With live hits it follows
  up on a random untried orthogonal neighbour of a live hit 3 times in 10; otherwise it fires at a random untried square.
  Modelled exactly (neighbour list keeps duplicates, as the Swift `flatMap` does).
- **Easy**: same with `followUpChance: 6` (ComputerPlayer.swift:33-34). Modelled exactly.
- **Normal**: `huntAndTarget` (ComputerPlayer.swift:35-36, 288-333). Follow-up: extend any line of two or more live
  hits from both ends, else any untried neighbour of a live hit. Search: only squares whose open run across or down
  fits the smallest ship afloat (`BotTargeting.run`, Opponents/BotTargeting.swift:132-146), preferring the pattern
  `(row + column) % smallest == 0`. Modelled exactly; the run test is computed as "some all-open window of that
  length covers the square", which is equivalent.
- **Hard**: `probabilityTarget(spreadsOpening: true)` with hunt-and-target as fallback (ComputerPlayer.swift:37-39,
  339-408). For every ship afloat, counts placements that avoid every marked square except live hits; with ships
  apart, squares diagonal to a live hit are water and a placement may not run alongside a live hit it doesn't
  include (lines 349-375). With live hits only placements through one count, weighted by hits covered squared
  (line 378). Fires at the best square (random tie-break, line 407), except that for the first 8 search shots
  (`openingSpreadShots`, line 337) it picks among the 5 best squares at random, weighted by score (lines 394-406).
  1-square ships are enumerated in both orientations, as in Swift. Modelled exactly.
- **Expert**: like Hard without the opening spread (`spreadsOpening: false`, ComputerPlayer.swift:40-42), and while
  searching it multiplies each square's score by `1 + 0.6 × habit` (lines 386-392). `habit` is the share of the
  player's past computer games with a ship on that square (GameModel.swift:3141-3162), learned only from the
  player's own fleet after each game ends (GameModel.swift:2840-2842). The main tables run Expert **with no history**
  (empty habits: identical to Hard minus the opening spread). The placement experiment adds
  **Expert after 20 games**: habits computed from 20 earlier fleets drawn from the same placement style, i.e. a
  player who always hides ships the same way.
- Not modelled: special weapons and defences of Advanced battles (`weaponPlan`, `defenses`), salvo rules, the Fog
  event rule, Repair Kit / Second Chance boosters, and automatic difficulty changes. None of them is part of a
  classic one-shot battle.
- Torpedo's own fleet placement for the computer and the Shuffle button (`Board.random(for:)`, Board.swift:46-63)
  places ships largest first, each uniformly among the spots still legal, and restarts on a dead end. The games
  here use uniform layouts instead, as requested; the heatmaps section shows how `Board.random` differs.
"""

OBSERVATIONS = """## Observations (written after the full run; numbers quoted from the tables above)

- **Torpedo's levels are well spread and its top levels match the textbook optimum.** Classic fleet means:
  Beginner 73.39, Easy 64.74, Normal 49.69, Hard 44.88, Expert 44.78 shots, against 45.05 for the published
  probability-density method. Hard's squared hit weighting gains about 0.2 shots on it; Hard's randomised opening
  (first 8 search shots) costs only about 0.1 shot against Expert with no history.
- **The port is faithful.** Against fleets from `Board.random`, the Python port and the real Swift code agree for
  all 10 level/fleet pairs within 1.4 standard errors.
- **Parity is worth about 7 shots** in the classic game (hunt/target 62.15 -> 55.13), and almost nothing in the
  10-ship game (63.18 -> 62.95), because a checkerboard is useless while any 1-square ship is afloat.
- **Random shooting needs 95.40 shots on average (median 97)** and never finished within 50 shots in 100,000 games.
- **Edges are safer only against a smart hunter, and only until it learns.** A fleet lying flat along the edges
  costs the probability-density hunter 51.32 shots instead of 45.10 (+6.2) and Torpedo Hard 50.66 instead of 44.76,
  but it is found *faster* by the hunt/target + parity hunter (53.43 vs 55.11), and Expert after 20 games against
  an edge-hugger needs only 46.80. Centre clusters are the worst choice (42.49 vs density; 36.27 vs a learned
  Expert), and packing all ships into one touching group also backfires (43.29 vs density, 47.22 vs parity).
- **The classic heatmap is milder than the naive count suggests.** Counting single-ship placements on an empty board
  gives a centre/corner ratio of 3.4; with all five ships placed together it is 2.67 (21.4% vs 8.0%).
- **Surprise: with the 10-ship no-touch fleet the map flips.** Edge squares are the hottest (best first shots
  C1, H1, A3, J3, A8, J8, C10, H10 at about 25.7%), the ring one square in from the edge is the coldest (14.7%),
  and the centre sits at 18.2%. A ship on the edge forbids fewer neighbouring squares, so many more legal layouts
  use the edges. The empty-board placement count, which puts the centre first (36 vs 16), points the wrong way here.
- **Torpedo's own random placer is not uniform.** `Board.random` places ships one by one: in the classic game it
  is a bit more centre-heavy than uniform (centre/corner 3.0 vs 2.67) and for the 10-ship fleet it is almost flat
  (1.06) instead of edge-heavy (0.77). The AI's results barely change (classic Hard 44.75 vs 44.88).
- Only 0.98% of uniformly random classic layouts have all five ships in one touching group.
- Run time: 11 to 13 minutes wall clock on a 14-core Mac with 12 worker processes (about 1 hour of CPU). The Swift
  cross-check harness (5,000 games per level, single thread) took about 9 minutes.
"""


if __name__ == "__main__":
    main()
