From 26de01b2e31d865cf876b4a4e0f362062c17fc4d Mon Sep 17 00:00:00 2001 From: Duncan Tourolle Date: Fri, 31 Jul 2026 22:28:06 +0200 Subject: [PATCH] study(VR-012): quality knee over a joint size x blur grid, three blur families Extends the VR-005 protocol -- hold out one mugshot per actor, degrade only the probe, match against a gallery held at native resolution, decide through the Platt calibration -- from one axis to two, over 1670 actors rather than 100. Joint rather than separable, because the interaction is the question: a 16 px face upscaled to 112 has already lost its high frequencies, so further blur costs it almost nothing, while the same blur at full resolution is expensive. Sweeping the axes independently would measure each with the other implicitly at its best and miss that entirely. Three blur families, compared at matched per-axis PSF spread rather than at equal raw parameter. Optical defocus is a uniform disc whose transfer function is a jinc with exact zeros, not a Gaussian that merely rolls off, and it is also how a face ends up large and useless -- the case a size filter cannot catch. Sweeping Gaussian alone, as the first version did, understates real lens blur by a factor of five in error rate. Every candidate measure is scored on every degraded crop and the candidates are ranked by how well each predicts the pipeline's actual decision, not by how smooth its synthetic ladder looks. Both a pooled and a within-cell AUC are reported: they answer different questions and the candidates rank differently under each. Runs through sae_embed throughout. Stages gains optional engine paths so the same study can drive a TRT build, which is what makes the full grid five minutes rather than four and a half hours. TRACES: VR-012, AR-028, AR-029 | SR-002 --- scripts/validation/min_face_size.py | 10 +- scripts/validation/quality_knee.py | 667 ++++++++++++++++++++++++++++ 2 files changed, 675 insertions(+), 2 deletions(-) create mode 100644 scripts/validation/quality_knee.py diff --git a/scripts/validation/min_face_size.py b/scripts/validation/min_face_size.py index 8c17538..74678b9 100644 --- a/scripts/validation/min_face_size.py +++ b/scripts/validation/min_face_size.py @@ -183,10 +183,16 @@ DEDUP_SIM = 1.0 - 1e-7 class Stages: """Thin holder so the rest of the script has one object to call.""" - def __init__(self, detector: str, arcface: str, conf: float, nms: float): + def __init__(self, detector: str, arcface: str, conf: float, nms: float, + detector_engine: str = "", arcface_engine: str = ""): + # The engine paths are only consulted by a TRT-backend build, where they + # are mandatory — that backend loads a pre-built .engine and will not + # fall back to reading the .onnx. An ORT build ignores them, so passing + # them unconditionally is safe and keeps one constructor for both. self.engine = sae_embed.FaceEmbedder( detector_model=detector, arcface_model=arcface, - conf=conf, nms=nms, max_side=0) + conf=conf, nms=nms, max_side=0, + detector_engine=detector_engine, arcface_engine=arcface_engine) def detect(self, img): return self.engine.detect(img) diff --git a/scripts/validation/quality_knee.py b/scripts/validation/quality_knee.py new file mode 100644 index 0000000..061cc14 --- /dev/null +++ b/scripts/validation/quality_knee.py @@ -0,0 +1,667 @@ +#!/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())