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())