#!/usr/bin/env python3
"""
quality_knee.py — VR-012: what does a blurred or small face cost in identification,
and which sharpness measure predicts it?

TRACES: VR-012, AR-028, AR-029

VR-005 located the size floor by degrading held-out gallery mugshots and watching
TPI/FPI fall. This does the same over a **joint size x blur grid**, and adds the
part that makes the result usable at inference.

Why a joint grid and not two sweeps
-----------------------------------
A 16 px face upscaled to 112 has already lost its high frequencies, so additional
blur costs it far less than it costs a 112 px one. Sweeping the axes separately
measures each in the presence of an implicit "other axis at its best" and misses
that interaction entirely — and the interaction is the whole question, because
AR-002 already gates on size and AR-029 proposes to discount on sharpness. If
identity loss turns out to be a function of the sharpness measure alone, then one
axis carries the information and discounting on both double-counts. If a
small-but-sharp and a large-but-blurred probe at equal measure lose different
amounts, the axes are genuinely separate and both belong.

Why sigma is not the answer
---------------------------
Sigma is a lab variable. At inference nothing knows how blurred a face is, so a
knee expressed in sigma cannot be acted on. What AR-028/AR-030 can consume is
    measure value -> expected identity reliability
so the controlled degradation exists to *select and calibrate the measure*, and
the measure is what ships. Every candidate is therefore scored on every degraded
crop, and the candidates are ranked by how well each predicts the identification
outcome (AUC over probe-cell records), not by how smooth its ladder looks.

Protocol (VR-005's, extended)
-----------------------------
1. Every gallery actor with at least `--min-images` mugshots. At the default 3,
   holding one out still leaves two references per actor.
2. Hold out ONE image per actor as the probe; the rest stay in the gallery at
   native resolution. Only the probe degrades — reference mugshots are clean and
   the face coming out of the video is not, which is the production case.
3. For each (size, sigma) cell: downscale the probe crop to size x size and back
   to 112 (the sampling loss), then Gaussian blur at sigma canonical px (the
   optical/motion loss). Resolution first, then blur, so sigma always means the
   same thing in the frame AR-029 measures in, whatever the cell's size.
4. Score all five AR-029 candidates on the degraded crop, through the C++
   binding.
5. Embed, match against the whole gallery, record TPI/FPI/unidentified.

Decision rule is the pipeline's: per-actor best-of-N cosine -> Platt sigmoid ->
accept if P > prob_threshold. Never a raw cosine (CLAUDE.md invariant, AR-024).

Everything runs through `sae_embed` — detection, the ArcFace warp, the embedder,
the sharpness measures and the calibration are all the shipped C++. Nothing here
re-implements a pipeline stage in numpy; the analysis on top of the recorded
numbers (AUC, knee location) is analysis and is numpy's job.

CAVEAT — FPI IS RELATIVE, NOT ABSOLUTE
--------------------------------------
False positives grow with the number of actors competing. Read FPI as a curve
across cells, not as a production rate. This runs the whole eligible gallery
rather than VR-005's 100-actor sample, so the understatement is much smaller,
but a production library is larger still.

Usage
-----
  python scripts/validation/quality_knee.py \
      --images images --gallery gallery_lvface.h5 \
      --arcface models/LVFace-B_Glint360K.onnx \
      --min-images 3 --out experiments/results/vr012_quality_knee
"""
from __future__ import annotations

import argparse
import csv
import json
import random
import sys
import time
from pathlib import Path

REPO = Path(__file__).resolve().parent.parent.parent
sys.path.insert(0, str(REPO / "scripts"))
sys.path.insert(0, str(Path(__file__).resolve().parent))

# min_face_size owns the shared scaffolding — actor discovery, the sae_embed
# locator, the Stages wrapper, the calibration-through-the-binding and the house
# plot palette. Importing it keeps one copy of each; a second copy of the
# calibration path in particular is what AR-024 exists to prevent.
import min_face_size as vr005          # noqa: E402
from min_face_size import (            # noqa: E402
    DEDUP_SIM, INTERP, Stages, calibrate_gallery, discover_actors,
    gallery_keys, normalise_name, probability, err,
    INK, MUTED, GRID, SURFACE, BLUE, GREEN, RED, AMBER,
)

import cv2        # noqa: E402
import numpy as np  # noqa: E402
import sae_embed   # noqa: E402

# The five AR-029 candidates, in quality.hpp's order. Names match the binding's
# attributes so the CSV columns and the C++ fields cannot drift apart.
MEASURES = ["var_laplacian", "norm_var_laplacian", "tenengrad",
            "hf_energy_ratio", "dir_min_tenengrad"]


# ── Degradation ───────────────────────────────────────────────────────────────

def disc_kernel(radius: float) -> np.ndarray:
    """The circle-of-confusion PSF of a defocused lens.

    Optical defocus is **not** Gaussian, and the difference is not cosmetic. A
    lens out of focus spreads a point into a uniform disc, whose transfer
    function is a jinc — `2·J1(x)/x` — which crosses zero and goes negative.
    Defocus therefore reverses contrast at particular spatial frequencies and
    can leave *more* energy in some high bands than a Gaussian of the same
    nominal width. A Gaussian MTF is strictly positive and monotonically
    decreasing and does neither.

    That matters here beyond realism: defocus is how a face ends up **large and
    useless**. A focus pull, a shallow depth of field, an actor stepping off the
    focal plane — all leave a big, confidently-detected face carrying no usable
    detail, and all sail straight through a size gate. Gaussian blur was the one
    family that mostly co-occurs with small faces, which is precisely why
    sharpness looked redundant against AR-002 on the first grid.

    The disc is supersampled 8x before downsampling so its edge is
    anti-aliased; a hard-edged binary disc at small radii is a poor circle and
    its spectrum carries the staircase, not the optics.
    """
    ss = 8
    n = int(np.ceil(radius)) * 2 + 1
    hi = np.zeros((n * ss, n * ss), np.float32)
    c = (n * ss - 1) / 2.0
    y, x = np.ogrid[:n * ss, :n * ss]
    hi[((x - c) ** 2 + (y - c) ** 2) <= (radius * ss) ** 2] = 1.0
    k = hi.reshape(n, ss, n, ss).mean(axis=(1, 3))
    s = k.sum()
    return (k / s) if s > 0 else np.ones((1, 1), np.float32)


def motion_kernel(length: int, angle_deg: float) -> np.ndarray:
    """Linear motion blur — a camera pan or a moving subject.

    Directional by construction: it destroys detail along one axis and leaves
    the perpendicular axis untouched. That is the property that separates the
    AR-029 candidates, since a measure normalising by total energy divides out
    the loss and reads a heavy smear as mild (see tests/test_quality.cpp).
    """
    k = np.zeros((length, length), np.float32)
    k[length // 2, :] = 1.0
    m = cv2.getRotationMatrix2D(((length - 1) / 2.0, (length - 1) / 2.0),
                                angle_deg, 1.0)
    k = cv2.warpAffine(k, m, (length, length))
    s = k.sum()
    return (k / s) if s > 0 else np.ones((1, 1), np.float32)


def degrade(crop: np.ndarray, size: int, level: float, kind: str,
            down: int, up: int, angle: float = 0.0) -> np.ndarray:
    """Resolution loss, then blur of the requested family.

    Order matters and this one is deliberate. Sampling happens in the source
    frame, so the downscale/upscale pair models a face that was `size` px when
    detected. The blur is then applied in the canonical frame, so `level` means
    the same number of canonical pixels in every cell of the grid — which is what
    lets the two axes be read independently. Blurring first would make the
    effective width depend on the cell's size, and the grid would no longer be
    factorial.

    `level` is the family's natural parameter: Gaussian sigma, disc radius, or
    motion length in canonical px. They are NOT equivalent at equal numbers —
    matching families by parameter would compare different amounts of damage, so
    the analysis matches them on measured effect instead.
    """
    out = crop
    if size != 112:
        small = cv2.resize(out, (size, size), interpolation=down)
        out = cv2.resize(small, (112, 112), interpolation=up)
    if level > 0:
        if kind == "gaussian":
            out = cv2.GaussianBlur(out, (0, 0), level, level)
        elif kind == "disc":
            out = cv2.filter2D(out, -1, disc_kernel(level))
        elif kind == "motion":
            out = cv2.filter2D(out, -1, motion_kernel(int(round(level)), angle))
        else:
            raise ValueError(f"unknown blur kind: {kind}")
    return out


def score_sharpness(crop: np.ndarray) -> dict:
    """All five candidates, from the shipped C++ (quality.hpp)."""
    s = sae_embed.assess_sharpness(np.ascontiguousarray(crop))
    d = {m: float(getattr(s, m)) for m in MEASURES}
    d["ok"] = bool(s.ok)
    return d


# ── Analysis ──────────────────────────────────────────────────────────────────

def auc(scores: np.ndarray, positive: np.ndarray) -> float:
    """Area under the ROC for `scores` predicting `positive`, by the rank
    (Mann-Whitney U) identity. 0.5 is chance; 1.0 is a measure that orders every
    correctly-identified probe above every failure.

    This is the ranking criterion for AR-029. A measure earns the job by
    predicting *the decision the pipeline makes*, not by having a tidy response
    to synthetic blur — a candidate can be beautifully monotone in sigma and
    still be a poor guide to whether this particular face will be recognised.
    """
    pos = scores[positive]
    neg = scores[~positive]
    if pos.size == 0 or neg.size == 0:
        return float("nan")
    order = np.argsort(np.concatenate([pos, neg]), kind="mergesort")
    ranks = np.empty(order.size, dtype=np.float64)
    ranks[order] = np.arange(1, order.size + 1)
    # Average ranks over ties, or a measure with many equal values is scored
    # arbitrarily by input order.
    vals = np.concatenate([pos, neg])
    sv = vals[order]
    i = 0
    while i < sv.size:
        j = i
        while j + 1 < sv.size and sv[j + 1] == sv[i]:
            j += 1
        if j > i:
            ranks[order[i:j + 1]] = ranks[order[i:j + 1]].mean()
        i = j + 1
    r_pos = ranks[:pos.size].sum()
    return float((r_pos - pos.size * (pos.size + 1) / 2) / (pos.size * neg.size))


def knee_from_measure(records: list[dict], measure: str, retention: float,
                      n_bins: int = 20) -> dict:
    """Where on `measure`'s own scale does identification start to fall apart?

    Bins the probe-cell records by measure value and reports the TPI rate in
    each. The threshold is the lowest bin edge whose bin and every bin above it
    retain `retention` of the undegraded control's TPI rate — a stated rule, so
    changing the answer means changing the rule rather than picking a number.
    """
    vals = np.array([r[measure] for r in records], dtype=np.float64)
    tpi = np.array([r["outcome"] == "TPI" for r in records])
    control = np.array([r["size_px"] == 112 and r["sigma"] == 0.0 for r in records])
    if control.sum() == 0:
        return {}
    floor = retention * float(tpi[control].mean())

    # Quantile edges: the measures have wildly different scales and heavy tails,
    # so equal-width bins would put almost everything in one bucket.
    edges = np.unique(np.quantile(vals, np.linspace(0, 1, n_bins + 1)))
    if edges.size < 3:
        return {}
    idx = np.clip(np.digitize(vals, edges[1:-1]), 0, edges.size - 2)

    bins = []
    for b in range(edges.size - 1):
        m = idx == b
        if m.sum() == 0:
            continue
        bins.append({"lo": float(edges[b]), "hi": float(edges[b + 1]),
                     "n": int(m.sum()), "tpi_rate": float(tpi[m].mean()),
                     "fpi_rate": float(np.mean([r["outcome"] == "FPI"
                                                for r, k in zip(records, m) if k]))})
    # Walk down from the top; the threshold is where retention first breaks.
    thr = None
    for b in reversed(bins):
        if b["tpi_rate"] < floor:
            thr = b["hi"]
            break
    return {"measure": measure, "control_tpi": float(tpi[control].mean()),
            "tpi_floor": floor, "threshold": thr, "bins": bins}


# ── Plot ──────────────────────────────────────────────────────────────────────

def write_plots(cells: list[dict], records: list[dict], ranking: list[dict],
                out_png: Path, meta: dict) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    plt.rcParams.update({
        "figure.facecolor": SURFACE, "axes.facecolor": SURFACE,
        "savefig.facecolor": SURFACE, "text.color": INK,
        "axes.edgecolor": MUTED, "axes.labelcolor": INK,
        "xtick.color": MUTED, "ytick.color": MUTED,
        "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.8,
        "axes.spines.top": False, "axes.spines.right": False,
    })

    sizes = sorted({c["size_px"] for c in cells})
    sigmas = sorted({c["sigma"] for c in cells})
    fig, axes = plt.subplots(1, 3, figsize=(17, 5.4))

    # (a) the joint grid as TPI heat map
    grid = np.full((len(sigmas), len(sizes)), np.nan)
    for c in cells:
        grid[sigmas.index(c["sigma"]), sizes.index(c["size_px"])] = 100 * c["tpi_rate"]
    im = axes[0].imshow(grid, origin="lower", aspect="auto", cmap="viridis",
                        vmin=0, vmax=100)
    axes[0].set_xticks(range(len(sizes)), [str(s) for s in sizes])
    axes[0].set_yticks(range(len(sigmas)), [f"{s:g}" for s in sigmas])
    axes[0].set_xlabel("probe size before upscaling (px)")
    axes[0].set_ylabel("Gaussian sigma (canonical px)")
    axes[0].set_title("TPI % over the joint grid", fontsize=11, loc="left")
    axes[0].grid(False)
    fig.colorbar(im, ax=axes[0], fraction=0.046)

    # (b) TPI against the winning measure — the curve a discount is built from
    best = ranking[0]["measure"]
    vals = np.array([r[best] for r in records])
    tpi = np.array([r["outcome"] == "TPI" for r in records])
    edges = np.unique(np.quantile(vals, np.linspace(0, 1, 21)))
    centres, rates = [], []
    for i in range(edges.size - 1):
        m = (vals >= edges[i]) & (vals <= edges[i + 1])
        if m.sum() > 20:
            centres.append(0.5 * (edges[i] + edges[i + 1]))
            rates.append(100 * tpi[m].mean())
    axes[1].plot(centres, rates, "-o", color=GREEN, lw=2)
    axes[1].set_xscale("log")
    axes[1].set_xlabel(f"{best} (log scale)")
    axes[1].set_ylabel("TPI %")
    axes[1].set_title(f"identification vs the measure\nbest predictor: {best} "
                      f"(AUC {ranking[0]['auc']:.3f})", fontsize=11, loc="left")

    # (c) how well each candidate predicts the decision
    names = [r["measure"] for r in ranking]
    aucs = [r["auc"] for r in ranking]
    axes[2].barh(range(len(names)), aucs, color=BLUE)
    axes[2].axvline(0.5, color=RED, lw=1.4, ls="--")
    axes[2].set_yticks(range(len(names)), names, fontsize=9)
    axes[2].set_xlim(0.4, 1.0)
    axes[2].set_xlabel("AUC — predicts correct identification")
    axes[2].set_title("AR-029 candidate ranking", fontsize=11, loc="left")
    axes[2].invert_yaxis()

    fig.suptitle(f"VR-012 — quality knee, {meta['model']}, {meta['n_actors']} actors, "
                 f"{meta['n_probes']} probes x {len(cells)} cells",
                 fontsize=12, x=0.01, ha="left")
    fig.tight_layout(rect=(0, 0.02, 1, 0.97))
    out_png.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(out_png, dpi=150)
    plt.close(fig)


# ── Main ──────────────────────────────────────────────────────────────────────

def main() -> int:
    p = argparse.ArgumentParser(
        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("--images", default=str(REPO / "images"))
    p.add_argument("--gallery", default=str(REPO / "gallery_lvface.h5"))
    p.add_argument("--out", default=str(REPO / "experiments/results/vr012_quality_knee"))
    p.add_argument("--actors", type=int, default=0,
                   help="cap the actor pool (0 = every eligible actor, the default: "
                        "FPI is gallery-size dependent and the whole gallery is the "
                        "least understated estimate available)")
    p.add_argument("--min-images", type=int, default=3,
                   help="minimum mugshots to be eligible (default 3, so holding one "
                        "out still leaves two references)")
    p.add_argument("--seed", type=int, default=0)
    p.add_argument("--sizes", default="16,24,32,48,64,112",
                   help="probe sizes before upscaling; 112 is undegraded")
    p.add_argument("--sigmas", default="0,0.5,1,1.5,2,3",
                   help="blur level in canonical px; 0 is unblurred. Meaning "
                        "depends on --blur-kind: Gaussian sigma, disc radius, "
                        "or motion length")
    p.add_argument("--blur-kind", default="gaussian",
                   choices=["gaussian", "disc", "motion"],
                   help="blur family. gaussian is a soft-focus stand-in; disc "
                        "is the circle-of-confusion PSF of real optical "
                        "defocus (non-Gaussian, jinc MTF with zero crossings); "
                        "motion is a linear smear. The last two are how a face "
                        "ends up large and useless, which a size gate cannot "
                        "catch")
    p.add_argument("--motion-angle", type=float, default=0.0,
                   help="motion blur direction in degrees (--blur-kind motion)")
    p.add_argument("--keep-duplicates", action="store_true")

    p.add_argument("--models-dir", default=str(REPO / "models"))
    p.add_argument("--arcface", default=None)
    p.add_argument("--detector", default=None)
    p.add_argument("--conf", type=float, default=0.5)
    p.add_argument("--nms", type=float, default=0.4)
    p.add_argument("--max-side", type=int, default=500)
    # Required by a TRT-backend build, ignored by an ORT one. A TensorRT fp16
    # run is a different realisation of the embedder — VR-005 measured ~0.85
    # cosine agreement with the fp32 ONNX path on LVFace-B, with separation
    # essentially intact — so a knee located here belongs to the fp16 space.
    # The study stays internally consistent because gallery and probes are both
    # embedded in this one session.
    p.add_argument("--detector-engine", default="",
                   help="pre-built SCRFD .engine (TRT builds only)")
    p.add_argument("--arcface-engine", default="",
                   help="pre-built ArcFace .engine (TRT builds only)")

    p.add_argument("--prob-threshold", type=float, default=0.754)
    p.add_argument("--match-prior", type=float, default=0.5)
    p.add_argument("--tpi-retention", type=float, default=0.95)
    p.add_argument("--down-interp", default="area", choices=sorted(INTERP))
    p.add_argument("--up-interp", default="linear", choices=sorted(INTERP))
    args = p.parse_args()

    models_dir = Path(args.models_dir)
    arcface = Path(args.arcface) if args.arcface else models_dir / "LVFace-B_Glint360K.onnx"
    detector = Path(args.detector) if args.detector else models_dir / "scrfd_500m_bnkps.onnx"
    for path, what in ((arcface, "embedder"), (detector, "detector")):
        if not path.is_file():
            return err(f"{what} model not found: {path}")

    images_root = Path(args.images)
    if not images_root.is_dir():
        return err(f"image cache not found: {images_root}")

    sizes = sorted({int(s) for s in args.sizes.split(",") if s.strip()})
    sigmas = sorted({float(s) for s in args.sigmas.split(",") if s.strip()})
    cv2.setRNGSeed(args.seed)

    # ── actor pool ────────────────────────────────────────────────────────────
    pool = discover_actors(images_root)
    print(f"[select] {len(pool)} actor dirs under {images_root}", file=sys.stderr)
    if args.gallery and Path(args.gallery).is_file():
        ids, names = gallery_keys(Path(args.gallery))
        pool = [a for a in pool
                if (a["jellyfin_id"] and a["jellyfin_id"] in ids)
                or normalise_name(a["name"]) in names]
        print(f"[select] {len(pool)} are in {args.gallery}", file=sys.stderr)

    eligible = [a for a in pool if len(a["images"]) >= args.min_images]
    print(f"[select] {len(eligible)} have >= {args.min_images} mugshots",
          file=sys.stderr)
    if len(eligible) < 2:
        return err(f"need at least 2 eligible actors; found {len(eligible)}")

    rng = random.Random(args.seed)
    selected = (sorted(rng.sample(eligible, min(args.actors, len(eligible))),
                       key=lambda a: a["dir"].name)
                if args.actors else eligible)

    # ── detect + align every mugshot once ─────────────────────────────────────
    stages = Stages(str(detector), str(arcface), args.conf, args.nms,
                    args.detector_engine, args.arcface_engine)
    print(f"[models] detector={detector.name} embedder={arcface.name} "
          f"batch={stages.engine.max_batch}", file=sys.stderr)

    t0 = time.time()
    crops, rows, actors = [], [], []
    n_nodetect = 0
    for a in selected:
        actor_crops, actor_paths = [], []
        for img_path in a["images"]:
            img = cv2.imread(str(img_path))
            if img is None:
                n_nodetect += 1
                continue
            if args.max_side > 0 and max(img.shape[:2]) > args.max_side:
                s = args.max_side / max(img.shape[:2])
                img = cv2.resize(img, None, fx=s, fy=s, interpolation=cv2.INTER_AREA)
            faces = stages.detect(img)
            if not faces:
                enhanced = stages.enhance(img)
                faces = stages.detect(enhanced)
                if faces:
                    img = enhanced
            if not faces:
                n_nodetect += 1
                continue
            best = max(faces, key=lambda f: f.confidence)
            crop = stages.align(img, best.landmarks)
            if crop is None:
                n_nodetect += 1
                continue
            actor_crops.append(crop)
            actor_paths.append(img_path)
        if len(actor_crops) < 2:
            continue
        ai = len(actors)
        actors.append({"name": a["name"], "jellyfin_id": a["jellyfin_id"],
                       "dir": a["dir"].name, "n_images": len(actor_crops)})
        for crop, img_path in zip(actor_crops, actor_paths):
            rows.append({"actor_idx": ai, "image": str(img_path)})
            crops.append(crop)
        if len(actors) % 200 == 0:
            print(f"  [align] {len(actors)}/{len(selected)} actors, "
                  f"{len(crops)} crops", file=sys.stderr)

    if len(actors) < 2:
        return err(f"only {len(actors)} actors survived detection/alignment")
    print(f"[align] {len(actors)} actors, {len(crops)} crops, {n_nodetect} skipped "
          f"in {time.time() - t0:.1f}s", file=sys.stderr)

    actor_of = np.array([r["actor_idx"] for r in rows], dtype=int)

    t0 = time.time()
    native = stages.embed(crops)
    print(f"[embed] {len(crops)} native crops in {time.time() - t0:.1f}s",
          file=sys.stderr)

    # ── drop duplicate mugshots ───────────────────────────────────────────────
    if not args.keep_duplicates:
        keep = np.ones(len(rows), bool)
        for ai in range(len(actors)):
            kept: list[int] = []
            for i in np.nonzero(actor_of == ai)[0]:
                if any(float(native[i] @ native[k]) > DEDUP_SIM for k in kept):
                    keep[i] = False
                else:
                    kept.append(int(i))
        n_dup = int((~keep).sum())
        counts = np.bincount(actor_of[keep], minlength=len(actors))
        drop_actor = counts < 2
        keep &= ~drop_actor[actor_of]
        remap = np.full(len(actors), -1, dtype=int)
        remap[~drop_actor] = np.arange(int((~drop_actor).sum()))
        actors = [a for a, d in zip(actors, drop_actor) if not d]
        rows = [r for r, k in zip(rows, keep) if k]
        crops = [c for c, k in zip(crops, keep) if k]
        native = native[keep]
        actor_of = remap[actor_of[keep]]
        print(f"[dedup] dropped {n_dup} duplicates and {int(drop_actor.sum())} "
              f"actors; {len(actors)} actors, {len(rows)} images remain",
              file=sys.stderr)

    # ── hold out one probe per actor ──────────────────────────────────────────
    is_probe = np.zeros(len(rows), bool)
    for ai in range(len(actors)):
        idx = np.nonzero(actor_of == ai)[0]
        r = random.Random(f"{args.seed}:{actors[ai]['dir']}")
        is_probe[r.choice(list(idx))] = True
    probe_rows = np.nonzero(is_probe)[0]
    gal_rows = np.nonzero(~is_probe)[0]
    print(f"[holdout] {len(probe_rows)} probes, {len(gal_rows)} gallery embeddings",
          file=sys.stderr)

    gal_emb = native[gal_rows]
    gal_actor = actor_of[gal_rows]
    probe_actor = actor_of[probe_rows]
    actor_cols = [np.nonzero(gal_actor == ai)[0] for ai in range(len(actors))]
    if not all(len(c) for c in actor_cols):
        return err("an actor has no gallery references left; raise --min-images")

    cal = calibrate_gallery(gal_emb, gal_actor)
    if not cal["valid"]:
        return err("calibration could not be fitted; this study will not fall back "
                   "to a raw cosine threshold (CLAUDE.md invariant)")
    log_prior_odds = float(np.log(args.match_prior / (1.0 - args.match_prior)))

    # ── the grid ──────────────────────────────────────────────────────────────
    down, up = INTERP[args.down_interp], INTERP[args.up_interp]
    probe_crops = [crops[i] for i in probe_rows]
    cells, records = [], []
    n = len(probe_rows)

    for size in sizes:
        for sigma in sigmas:
            t0 = time.time()
            degraded = [degrade(c, size, sigma, args.blur_kind, down, up,
                                args.motion_angle) for c in probe_crops]
            sharp = [score_sharpness(d) for d in degraded]
            q = stages.embed(degraded)

            sims = q @ gal_emb.T
            best_per_actor = np.stack([sims[:, c].max(axis=1) for c in actor_cols],
                                      axis=1)
            best_actor = best_per_actor.argmax(axis=1)
            best_sim = best_per_actor.max(axis=1)
            p_match = np.asarray(probability(best_sim, cal["a"], cal["b"],
                                             log_prior_odds))
            accept = p_match > args.prob_threshold
            correct = best_actor == probe_actor
            tpi = int(np.sum(accept & correct))
            fpi = int(np.sum(accept & ~correct))
            unid = int(np.sum(~accept))

            cell = {"size_px": size, "sigma": sigma, "blur_kind": args.blur_kind,
                    "n_probes": n,
                    "tpi": tpi, "fpi": fpi, "unidentified": unid,
                    "tpi_rate": tpi / n, "fpi_rate": fpi / n,
                    "unidentified_rate": unid / n,
                    "rank1_rate": float(np.mean(correct)),
                    "mean_p_match": float(np.mean(p_match))}
            for m in MEASURES:
                cell[f"mean_{m}"] = float(np.mean([s[m] for s in sharp]))
            cells.append(cell)

            for j in range(n):
                rec = {"size_px": size, "sigma": sigma,
                       "blur_kind": args.blur_kind,
                       "probe_image": rows[probe_rows[j]]["image"],
                       "p_match": float(p_match[j]),
                       "outcome": ("TPI" if accept[j] and correct[j]
                                   else "FPI" if accept[j] else "unidentified")}
                rec.update({m: sharp[j][m] for m in MEASURES})
                records.append(rec)

            print(f"[grid] {size:3d}px {args.blur_kind[:4]} {sigma:<4g}  TPI {100*tpi/n:5.1f}%  "
                  f"FPI {100*fpi/n:5.1f}%  unid {100*unid/n:5.1f}%  "
                  f"rank1 {100*np.mean(correct):5.1f}%  [{time.time()-t0:.1f}s]",
                  file=sys.stderr)

    # ── rank the candidates, then locate the knee on the winner ───────────────
    is_tpi = np.array([r["outcome"] == "TPI" for r in records])
    ranking = sorted(
        ({"measure": m,
          "auc": auc(np.array([r[m] for r in records], dtype=np.float64), is_tpi)}
         for m in MEASURES),
        key=lambda d: -d["auc"])
    knees = [knee_from_measure(records, r["measure"], args.tpi_retention)
             for r in ranking]

    # ── outputs ───────────────────────────────────────────────────────────────
    out = Path(args.out)
    out.parent.mkdir(parents=True, exist_ok=True)
    csv_path = out.with_name(out.name + ".csv")
    with open(csv_path, "w", newline="") as f:
        w = csv.DictWriter(f, fieldnames=list(cells[0].keys()))
        w.writeheader()
        w.writerows(cells)
    rec_path = out.with_name(out.name + ".records.csv")
    with open(rec_path, "w", newline="") as f:
        w = csv.DictWriter(f, fieldnames=list(records[0].keys()))
        w.writeheader()
        w.writerows(records)

    meta = {
        "requirement": "VR-012",
        "model": arcface.stem, "detector": detector.stem,
        "n_actors": len(actors), "n_probes": len(probe_rows),
        "n_gallery_embeddings": len(gal_rows),
        "min_images": args.min_images, "seed": args.seed,
        "sizes": sizes, "sigmas": sigmas, "blur_kind": args.blur_kind,
        "motion_angle": args.motion_angle,
        "prob_threshold": args.prob_threshold, "match_prior": args.match_prior,
        "calibration": cal,
        "measure_ranking": ranking,
        "knees": knees,
        "sharpness_window": list(sae_embed.sharpness_window()),
        "caveat": (f"FPI grows with gallery size; this ran against {len(actors)} "
                   f"actors and still understates a production library."),
        "grid": cells,
    }
    json_path = out.with_name(out.name + ".json")
    json_path.write_text(json.dumps(meta, indent=2) + "\n")
    png_path = out.with_name(out.name + ".png")
    write_plots(cells, records, ranking, png_path, meta)

    # ── stdout report ─────────────────────────────────────────────────────────
    print(f"\nVR-012 — quality knee, {arcface.stem}")
    print(f"{len(actors)} actors, {len(probe_rows)} probes x {len(cells)} cells\n")
    print(f"{'size':>5} {'sigma':>6} {'TPI':>8} {'FPI':>8} {'unid':>8} {'rank1':>8}")
    for c in cells:
        print(f"{c['size_px']:>5} {c['sigma']:>6g} {100*c['tpi_rate']:>7.1f}% "
              f"{100*c['fpi_rate']:>7.1f}% {100*c['unidentified_rate']:>7.1f}% "
              f"{100*c['rank1_rate']:>7.1f}%")
    print("\nAR-029 candidate ranking — AUC for predicting correct identification:")
    for r in ranking:
        print(f"  {r['measure']:>20}  {r['auc']:.4f}")
    print(f"\n[out] {csv_path}\n[out] {rec_path}\n[out] {json_path}\n[out] {png_path}")
    print(f"\n{meta['caveat']}")
    return 0


if __name__ == "__main__":
    sys.exit(main())