rpx_benchmark.data.esd_scoring

  1"""ESD scoring methods: turn the 27-feature table into a single difficulty
  2score per (scene, phase) and an Easy/Medium/Hard tertile assignment.
  3
  4Five scoring methods are exposed, all operating on the **percentile-normalised**
  5feature matrix (per-feature ``rank/N``, ∈ ``[0,1]``):
  6
  7- ``mean_pn``       — uniform weighted mean (paper-stated method, primary default)
  8- ``median_pn``     — robust to outlier features
  9- ``max_pn``        — weakest-link: hardest single feature defines difficulty
 10- ``pca_pc1``       — first principal component score
 11- ``kmeans3``       — k-means clusters tier-ordered by centroid magnitude
 12
 13Each method is a pure function that takes the percentile-normalised matrix
 14(plus optional weights for ``mean_pn``) and returns a ``(continuous_score,
 15categorical_tier)`` pair.
 16
 17Two weighting tiers are supported (paper §3.2):
 18
 19- ``uniform_v1`` — ``wᵢ = 1/F``. Always reproducible from data alone. The
 20  test-of-time baseline.
 21- ``mi_v1``     — weights derived from mutual information between features and
 22  model failure rates. Loaded from a JSON artifact (see :func:`load_mi_weights`).
 23  The headline "official" tier when calibration models are available.
 24
 25Heavy deps (``scikit-learn``) are lazy-imported inside each function so the
 26module imports cleanly on a minimal install — the ``mean_pn``, ``median_pn``,
 27``max_pn`` paths work without sklearn.
 28"""
 29
 30from __future__ import annotations
 31
 32import hashlib
 33import json
 34from dataclasses import dataclass, field
 35from pathlib import Path
 36from typing import Dict, List, Optional, Tuple
 37
 38import numpy as np
 39
 40from ..exceptions import ConfigError, DatasetError
 41from ..logging_utils import get_logger
 42
 43log = get_logger(__name__)
 44
 45
 46# Method names — kept as a tuple so callers can iterate stably.
 47SCORING_METHODS: Tuple[str, ...] = (
 48    "mean_pn",
 49    "median_pn",
 50    "max_pn",
 51    "pca_pc1",
 52    "kmeans3",
 53)
 54
 55# The default primary method. Documented in the paper as the v1 official score.
 56DEFAULT_PRIMARY: str = "mean_pn"
 57
 58# Tertile labels — kept as a constant so the order matches every output file.
 59TERTILE_LABELS: Tuple[str, str, str] = ("easy", "medium", "hard")
 60
 61# Confidence labels for per-row confidence based on weight-perturbation flip rate.
 62CONFIDENCE_LABELS: Tuple[str, str, str] = ("high", "medium", "low")
 63
 64# Feature → modality category mapping. Keys ordered by the paper appendix
 65# enumeration; values list the canonical FEATURE_NAMES that belong to each
 66# category. Used to compute the per-category sub-score vector that surfaces
 67# the multi-faceted nature of difficulty (data story §S4).
 68FEATURE_CATEGORIES: Dict[str, Tuple[str, ...]] = {
 69    "annotation_effort": ("iter_mean", "iter_max"),
 70    "scene_complexity": ("obj_mean", "obj_std", "obj_consist"),
 71    "occlusion": ("occ_mean", "occ_p90", "occ_heavy"),
 72    "depth_quality": ("depth_invalid", "depth_invalid_mask", "depth_std", "depth_std_mask"),
 73    "photometric_conflict": ("specular", "dark"),
 74    "image_quality": ("rgb_blur", "rgb_texture"),
 75    "object_size": ("mask_area_mean", "mask_area_std"),
 76    "temporal_stability": ("area_cv", "area_drop", "vis_instability"),
 77    "camera_motion": ("trans_mean", "trans_p90", "rot_mean", "rot_p90", "jerk"),
 78    "fisheye_stereo": (
 79        "fisheye_dark",
 80        "fisheye_bright",
 81        "fisheye_sharpness",
 82        "fisheye_corr",
 83        "fisheye_texture",
 84    ),
 85}
 86
 87
 88# --------------------------------------------------------------------------- #
 89# Provenance
 90# --------------------------------------------------------------------------- #
 91
 92
 93@dataclass
 94class ScoringProvenance:
 95    """Metadata stamped on every output file so a 2030 researcher can verify
 96    reproducibility without ambiguity."""
 97
 98    input_sha256: str
 99    feature_set_version: str
100    scoring_version: str  # "uniform_v1" | "mi_v1"
101    primary_method: str  # one of SCORING_METHODS
102    weights: Dict[str, float]  # {feature_name: weight} for primary scoring
103    feature_names: List[str]  # the 27 feature names, in canonical order
104    n_entries: int
105    generated_at_utc: str
106    script_version: str = "unknown"  # caller fills with `git rev-parse HEAD`
107    extra: Dict[str, object] = field(default_factory=dict)
108
109    def to_dict(self) -> Dict[str, object]:
110        return {
111            "input_sha256": self.input_sha256,
112            "feature_set_version": self.feature_set_version,
113            "scoring_version": self.scoring_version,
114            "primary_method": self.primary_method,
115            "weights": self.weights,
116            "feature_names": self.feature_names,
117            "n_entries": self.n_entries,
118            "generated_at_utc": self.generated_at_utc,
119            "script_version": self.script_version,
120            "extra": self.extra,
121        }
122
123
124def sha256_of_file(path: Path) -> str:
125    """SHA-256 hex digest of a file's bytes — used to anchor provenance."""
126    h = hashlib.sha256()
127    with Path(path).open("rb") as f:
128        for chunk in iter(lambda: f.read(64 * 1024), b""):
129            h.update(chunk)
130    return h.hexdigest()
131
132
133# --------------------------------------------------------------------------- #
134# Weights
135# --------------------------------------------------------------------------- #
136
137
138def uniform_weights(feature_names: List[str]) -> Dict[str, float]:
139    """``wᵢ = 1/F`` — the always-reproducible baseline.
140
141    Note: this is the ``uniform_v1`` tier. The methodology's headline scoring
142    is ``effort_stratified_v1`` (see :func:`effort_stratified_weights`) which
143    actually honours the methodology's name; ``uniform_v1`` is kept as a
144    "no-prior" reference point for the sensitivity analysis.
145    """
146    if not feature_names:
147        raise ConfigError(
148            "uniform_weights requires a non-empty feature list",
149            hint="check that FEATURE_NAMES is populated.",
150        )
151    w = 1.0 / len(feature_names)
152    return {name: w for name in feature_names}
153
154
155# Features that anchor the Effort-Stratified Difficulty (ESD) score: mask
156# refinement iteration counts. The methodology is *named* after these — the
157# ``effort_stratified_v1`` weighting reflects that name in the actual score.
158EFFORT_FEATURES: Tuple[str, ...] = ("iter_mean", "iter_max")
159
160# Default convex-combination weight on the effort block. α=0.25 means
161# "effort gets a modest boost" without dominating: each iter feature
162# contributes 12.5% of the score (4× the 3% each other feature gets),
163# but the 25 perception features still collectively drive 75% of the score.
164DEFAULT_EFFORT_ALPHA: float = 0.25
165
166
167def effort_stratified_weights(
168    feature_names: List[str],
169    *,
170    alpha: float = DEFAULT_EFFORT_ALPHA,
171    effort_features: Tuple[str, ...] = EFFORT_FEATURES,
172) -> Dict[str, float]:
173    """Effort-anchored weights — the headline ``effort_stratified_v1`` scoring.
174
175    The methodology is named *Effort-Stratified Difficulty*; this weighting
176    actually puts effort at the centre. The convex combination
177
178        RPX-DS = α · effort_score + (1 − α) · perception_score
179
180    expands to per-feature weights:
181
182        iter_mean = iter_max = α / |effort_features|
183        other features       = (1 − α) / (F − |effort_features|)
184
185    With ``alpha = 0.25`` (default) and 31 features (2 effort + 29 other),
186    iter features each carry 12.5 % of the score (25 % combined); each of the
187    29 perception features carries ~2.6 % (75 % combined). At ``alpha = 0`` the
188    weighting reduces to uniform over the perception features only; at
189    ``alpha = 1`` the score depends only on annotation effort.
190    """
191    if not feature_names:
192        raise ConfigError(
193            "effort_stratified_weights requires a non-empty feature list",
194            hint="check that FEATURE_NAMES is populated.",
195        )
196    if not (0.0 <= alpha <= 1.0):
197        raise ConfigError(
198            f"alpha must be in [0, 1], got {alpha}",
199            hint="alpha controls the convex combination of effort vs perception.",
200        )
201
202    effort = [n for n in feature_names if n in effort_features]
203    other = [n for n in feature_names if n not in effort_features]
204    if not effort:
205        raise ConfigError(
206            f"no effort features found in feature list (looked for {effort_features})",
207            hint="check that FEATURE_NAMES includes iter_mean / iter_max.",
208        )
209    if not other:
210        # Degenerate edge case: every feature is an effort feature → fall back
211        # to uniform within the effort block so weights still sum to 1.
212        return uniform_weights(feature_names)
213
214    w_effort = alpha / len(effort)
215    w_other = (1.0 - alpha) / len(other)
216    weights = {n: (w_effort if n in effort_features else w_other) for n in feature_names}
217    # Sanity: weights sum to 1 (modulo floating-point).
218    total = sum(weights.values())
219    assert abs(total - 1.0) < 1e-9, f"weights sum to {total}, expected 1.0"
220    return weights
221
222
223def load_mi_weights(
224    path: Path,
225    feature_names: List[str],
226) -> Tuple[Dict[str, float], Dict[str, object]]:
227    """Load mutual-information-derived weights from a JSON artifact.
228
229    The artifact contract::
230
231        {
232          "weights": { "<feature_name>": <float>, ... },
233          "calibration_model_set": [ "model_id", ... ],
234          "metric_used": "MI",
235          "notes": "..."
236        }
237
238    Weights are renormalised to sum to 1.0 over the requested ``feature_names``.
239    Features missing from the artifact get weight 0.0; features in the artifact
240    but not in ``feature_names`` are silently dropped (forward-compat: future
241    feature additions don't break old MI files).
242    """
243    path = Path(path)
244    if not path.is_file():
245        raise DatasetError(
246            f"MI weights file not found: {path}",
247            hint="pass --weights-mi to point at the calibration artifact.",
248        )
249    payload = json.loads(path.read_text())
250    raw = payload.get("weights")
251    if not isinstance(raw, dict):
252        raise DatasetError(
253            f"MI weights file missing 'weights' dict: {path}",
254            hint="artifact must look like {'weights': {feature: float, ...}, ...}",
255        )
256
257    aligned = {name: float(raw.get(name, 0.0)) for name in feature_names}
258    total = sum(aligned.values())
259    if total <= 0:
260        raise DatasetError(
261            f"MI weights sum to {total} after alignment with feature list",
262            hint="check that weights file uses the same feature names as the build.",
263        )
264    aligned = {k: v / total for k, v in aligned.items()}
265
266    metadata = {k: v for k, v in payload.items() if k != "weights"}
267    return aligned, metadata
268
269
270# --------------------------------------------------------------------------- #
271# Percentile-normalisation
272# --------------------------------------------------------------------------- #
273
274
275def percentile_normalize(matrix: np.ndarray) -> np.ndarray:
276    """Per-column rank divided by N (matches paper §3.2 ``f̃ᵢ``).
277
278    ``matrix`` shape ``(N, F)``. Output same shape, values in ``[0, 1]``.
279    Constant columns (zero variance) are mapped to all-0.5 — informationless
280    but non-degenerate; downstream scoring handles them gracefully.
281
282    NaN values (e.g., fisheye features for scenes without fisheye data) are
283    ranked only among non-NaN entries, then set to 0.5 (midpoint) in the
284    output. This prevents missing-modality scenes from being systematically
285    biased toward low difficulty.
286    """
287    if matrix.ndim != 2:
288        raise ConfigError(
289            f"percentile_normalize expects 2D array, got shape {matrix.shape}",
290            hint="reshape to (n_rows, n_features) before calling.",
291        )
292    n_rows, n_cols = matrix.shape
293    out = np.empty_like(matrix, dtype=np.float64)
294    for j in range(n_cols):
295        col = matrix[:, j]
296        valid_mask = ~np.isnan(col)
297        valid = col[valid_mask]
298        n_valid = valid.size
299
300        if n_valid == 0 or np.unique(valid).size <= 1:
301            out[:, j] = 0.5  # all-NaN or constant → uninformative midpoint
302            continue
303
304        # Rank only the non-NaN entries.
305        order = np.argsort(valid, kind="mergesort")
306        ranks = np.empty(n_valid, dtype=np.float64)
307        ranks[order] = np.arange(1, n_valid + 1, dtype=np.float64)
308        # Average ranks for tied values.
309        for v in np.unique(valid):
310            idx = np.where(valid == v)[0]
311            if idx.size > 1:
312                ranks[idx] = ranks[idx].mean()
313
314        # Write ranked values back; NaN entries get 0.5 (neutral midpoint).
315        out[:, j] = 0.5
316        out[valid_mask, j] = ranks / n_valid
317    return out
318
319
320# --------------------------------------------------------------------------- #
321# Tertile cut
322# --------------------------------------------------------------------------- #
323
324
325def tertile_cut(scores: np.ndarray) -> np.ndarray:
326    """Bottom 33% / middle 33% / top 34% → 'easy' / 'medium' / 'hard'.
327
328    Uses rank-based percentiles so ties never collapse a tertile. If the score
329    has fewer than 3 unique values everything degenerates to 'medium' (caller
330    can decide what to do).
331    """
332    n = scores.size
333    if n == 0:
334        return np.array([], dtype=object)
335    if np.unique(scores).size < 3:
336        log.warning("tertile_cut: fewer than 3 unique scores; defaulting all to 'medium'")
337        return np.array([TERTILE_LABELS[1]] * n, dtype=object)
338
339    ranks = np.argsort(np.argsort(scores, kind="mergesort")) + 1
340    p = ranks / n
341    out = np.empty(n, dtype=object)
342    out[p <= 1.0 / 3.0] = TERTILE_LABELS[0]
343    out[(p > 1.0 / 3.0) & (p <= 2.0 / 3)] = TERTILE_LABELS[1]
344    out[p > 2.0 / 3.0] = TERTILE_LABELS[2]
345    return out
346
347
348def percentile_rank(scores: np.ndarray) -> np.ndarray:
349    """Per-row percentile rank in ``[0, 1]`` based on ``argsort``."""
350    n = scores.size
351    if n == 0:
352        return np.array([], dtype=np.float64)
353    ranks = np.argsort(np.argsort(scores, kind="mergesort")) + 1
354    return ranks.astype(np.float64) / n
355
356
357# --------------------------------------------------------------------------- #
358# Scoring methods (each returns continuous score; tertile_cut applied separately)
359# --------------------------------------------------------------------------- #
360
361
362def score_mean_pn(pn: np.ndarray, weights: np.ndarray) -> np.ndarray:
363    """Weighted mean of percentile-normalised features. Default primary method.
364
365    With ``weights = 1/F`` this is the uniform RPX-DS in paper §3.2. With
366    MI-derived weights it becomes ``mi_v1``.
367    """
368    if pn.shape[1] != weights.size:
369        raise ConfigError(
370            f"weights size {weights.size} does not match feature count {pn.shape[1]}",
371            hint="weights must align with the feature columns of pn.",
372        )
373    return pn @ weights
374
375
376def score_median_pn(pn: np.ndarray) -> np.ndarray:
377    """Per-row median over features. Robust to outlier features."""
378    return np.median(pn, axis=1)
379
380
381def score_max_pn(pn: np.ndarray) -> np.ndarray:
382    """Per-row max over features. Weakest-link: hardest single feature wins."""
383    return pn.max(axis=1)
384
385
386def score_pca_pc1(pn: np.ndarray) -> np.ndarray:
387    """First principal component score of the centred percentile matrix.
388
389    PC1 is sign-flipped if needed so that higher score correlates with higher
390    mean_pn — keeps the "Hard = high score" convention consistent across methods.
391    """
392    try:
393        from sklearn.decomposition import PCA  # noqa: PLC0415
394    except ImportError as e:
395        raise ConfigError(
396            "score_pca_pc1 requires scikit-learn",
397            hint="install with: pip install 'rpx-benchmark[analysis]'",
398        ) from e
399    pca = PCA(n_components=1)
400    pc1 = pca.fit_transform(pn - pn.mean(axis=0)).ravel()
401    if np.corrcoef(pc1, pn.mean(axis=1))[0, 1] < 0:
402        pc1 = -pc1
403    return pc1
404
405
406def score_kmeans3(pn: np.ndarray, *, seed: int = 0) -> np.ndarray:
407    """K-means cluster id (k=3) on the first 5 PCs of the percentile matrix.
408
409    Returned as a *continuous* score: each row's value is the centroid's
410    rank by mean_pn-magnitude (0 → easiest cluster, 2 → hardest). The
411    ``tertile_cut`` then collapses these onto ``easy/medium/hard``.
412    """
413    try:
414        from sklearn.cluster import KMeans  # noqa: PLC0415
415        from sklearn.decomposition import PCA  # noqa: PLC0415
416    except ImportError as e:
417        raise ConfigError(
418            "score_kmeans3 requires scikit-learn",
419            hint="install with: pip install 'rpx-benchmark[analysis]'",
420        ) from e
421    n_rows = pn.shape[0]
422    if n_rows < 3:
423        log.warning("score_kmeans3: <3 rows; degenerate clustering, returning zeros")
424        return np.zeros(n_rows, dtype=np.float64)
425
426    n_pcs = min(5, pn.shape[1], n_rows)
427    proj = PCA(n_components=n_pcs).fit_transform(pn - pn.mean(axis=0))
428    km = KMeans(n_clusters=3, n_init=10, random_state=seed).fit(proj)
429    labels = km.labels_
430
431    # Order clusters by their mean of mean_pn so cluster 0 = easiest.
432    mean_pn = pn.mean(axis=1)
433    cluster_means = {c: float(mean_pn[labels == c].mean()) for c in np.unique(labels)}
434    sorted_clusters = sorted(cluster_means, key=cluster_means.get)
435    rank_of = {c: i for i, c in enumerate(sorted_clusters)}
436    return np.asarray([rank_of[c] for c in labels], dtype=np.float64)
437
438
439# --------------------------------------------------------------------------- #
440# Locality Preserving Projections (LPP) — He & Niyogi, NeurIPS 2003
441# --------------------------------------------------------------------------- #
442
443
444def lpp_embedding(
445    X: np.ndarray,
446    *,
447    n_components: int = 2,
448    n_neighbors: int = 10,
449    heat_t: Optional[float] = None,
450) -> np.ndarray:
451    """Locality Preserving Projections — linear projection that preserves
452    local neighbourhood structure on the data manifold.
453
454    Implements the generalised-eigenvalue formulation from
455    He & Niyogi (NeurIPS 2003): solve ``X^T L X v = λ X^T D X v`` where
456    ``L = D − W`` is the graph Laplacian over a kNN affinity matrix with
457    heat-kernel weights.
458
459    Parameters
460    ----------
461    X            : ``(n, d)`` data matrix
462    n_components : output dimensionality
463    n_neighbors  : k for the kNN affinity graph (excluding self)
464    heat_t       : heat-kernel bandwidth; defaults to median squared
465                   inter-neighbour distance
466
467    Returns
468    -------
469    ``(n, n_components)`` linear embedding ``Y = X @ W``.
470
471    Notes
472    -----
473    LPP is the linear approximation of Laplacian Eigenmaps. Compared to
474    sklearn's ``SpectralEmbedding`` (the non-linear version), LPP gives a
475    deterministic projection that can be applied to held-out data — useful
476    for analyses that need the projection matrix itself, not just the
477    embedding of the training points.
478    """
479    try:
480        from scipy.linalg import eigh  # noqa: PLC0415
481        from sklearn.neighbors import NearestNeighbors  # noqa: PLC0415
482    except ImportError as e:
483        raise ConfigError(
484            "lpp_embedding requires scikit-learn and scipy",
485            hint="install with: pip install 'rpx-benchmark[analysis]'",
486        ) from e
487
488    if X.ndim != 2:
489        raise ConfigError(
490            f"lpp_embedding expects 2D array, got shape {X.shape}",
491            hint="reshape to (n_rows, n_features) before calling.",
492        )
493    n, d = X.shape
494    if n_neighbors >= n:
495        raise ConfigError(
496            f"n_neighbors ({n_neighbors}) must be < n_samples ({n})",
497            hint="reduce n_neighbors or supply a larger dataset.",
498        )
499
500    # Step 1: kNN graph (k+1 to drop self).
501    nn = NearestNeighbors(n_neighbors=n_neighbors + 1).fit(X)
502    distances, indices = nn.kneighbors(X)
503    distances = distances[:, 1:]  # drop self
504    indices = indices[:, 1:]
505
506    # Step 2: heat-kernel weights, symmetrised affinity matrix W.
507    if heat_t is None:
508        heat_t = float(np.median(distances**2))
509        if heat_t <= 0:
510            heat_t = 1.0  # degenerate fallback (e.g., duplicated points)
511
512    W = np.zeros((n, n), dtype=np.float64)
513    for i in range(n):
514        for j_idx, j in enumerate(indices[i]):
515            w = float(np.exp(-(distances[i, j_idx] ** 2) / heat_t))
516            # Symmetrise: take max so the affinity is consistent under
517            # asymmetric kNN.
518            W[i, j] = max(W[i, j], w)
519            W[j, i] = max(W[j, i], w)
520
521    # Step 3: Laplacian.
522    deg = W.sum(axis=1)
523    D = np.diag(deg)
524    L = D - W
525
526    # Step 4: solve generalised eigenvalue problem X^T L X v = λ X^T D X v.
527    A = X.T @ L @ X
528    B = X.T @ D @ X
529    # Tikhonov regularisation on B to handle rank-deficient cases.
530    B = B + 1e-9 * np.eye(d)
531    eigvals, eigvecs = eigh(A, B)
532
533    # Take the n_components smallest non-zero eigenvalues (they correspond
534    # to directions that minimise local distortion).
535    proj = eigvecs[:, :n_components]
536    return X @ proj
537
538
539# --------------------------------------------------------------------------- #
540# Driver — run all methods at once
541# --------------------------------------------------------------------------- #
542
543
544def all_methods(
545    pn: np.ndarray,
546    weights: np.ndarray,
547    *,
548    seed: int = 0,
549) -> Dict[str, np.ndarray]:
550    """Compute every method's continuous score on the same percentile matrix.
551
552    Returns ``{method_name: 1D score array}``. Methods that need scikit-learn
553    are skipped with a warning if the import fails — the caller can then ship
554    the available subset.
555    """
556    out: Dict[str, np.ndarray] = {}
557    out["mean_pn"] = score_mean_pn(pn, weights)
558    out["median_pn"] = score_median_pn(pn)
559    out["max_pn"] = score_max_pn(pn)
560    try:
561        out["pca_pc1"] = score_pca_pc1(pn)
562    except ConfigError as e:
563        log.warning("skipping pca_pc1: %s", e)
564    try:
565        out["kmeans3"] = score_kmeans3(pn, seed=seed)
566    except ConfigError as e:
567        log.warning("skipping kmeans3: %s", e)
568    return out
569
570
571# --------------------------------------------------------------------------- #
572# Structured-metric helpers (data story §S4 + §S5 + consensus)
573# --------------------------------------------------------------------------- #
574
575# Per-row tier-flip rate thresholds for the confidence label. A pair flipping
576# tiers in <10% of weight perturbations is "high" confidence; 10-50% medium;
577# >50% low.
578DEFAULT_CONFIDENCE_HIGH_MAX: float = 0.10
579DEFAULT_CONFIDENCE_MEDIUM_MAX: float = 0.50
580
581
582def compute_per_category_scores(
583    pn: np.ndarray,
584    feature_names: List[str],
585    categories: Optional[Dict[str, Tuple[str, ...]]] = None,
586) -> Dict[str, np.ndarray]:
587    """Per-row mean of percentile-normalised features within each category.
588
589    Returns ``{category_name: 1D float array of length N}``.
590    """
591    cats = categories if categories is not None else FEATURE_CATEGORIES
592    name_to_idx = {n: i for i, n in enumerate(feature_names)}
593    out: Dict[str, np.ndarray] = {}
594    for cat, feats in cats.items():
595        idx = [name_to_idx[f] for f in feats if f in name_to_idx]
596        if not idx:
597            log.warning("category %s has no features in the matrix; skipping", cat)
598            continue
599        out[cat] = pn[:, idx].mean(axis=1)
600    return out
601
602
603def compute_consensus_tier(
604    method_tiers: Dict[str, np.ndarray],
605    *,
606    exclude_methods: Optional[Tuple[str, ...]] = None,
607) -> np.ndarray:
608    """Majority-vote tier across multiple methods' tier assignments.
609
610    Each method's tier array must have the same length. ``exclude_methods``
611    drops degenerate methods (e.g. DBSCAN on data with no density modes)
612    from the vote. Ties are broken by the canonical TERTILE_LABELS order
613    (easy < medium < hard) — i.e. ties favour the easier label, the
614    conservative choice for downstream evaluation.
615    """
616    from collections import Counter
617
618    excluded = set(exclude_methods or ())
619    voting = {m: t for m, t in method_tiers.items() if m not in excluded}
620    if not voting:
621        raise ConfigError(
622            "compute_consensus_tier received no methods to vote on",
623            hint="check that exclude_methods does not exhaust all methods.",
624        )
625    n = next(iter(voting.values())).shape[0]
626    out = np.empty(n, dtype=object)
627    for i in range(n):
628        votes = Counter(voting[m][i] for m in voting)
629        max_count = max(votes.values())
630        winners = [t for t, c in votes.items() if c == max_count]
631        if len(winners) == 1:
632            out[i] = winners[0]
633        else:
634            for t in TERTILE_LABELS:
635                if t in winners:
636                    out[i] = t
637                    break
638            else:
639                out[i] = winners[0]  # fallback (e.g. degenerate label)
640    return out
641
642
643def aggregate_to_scene_splits(
644    scene_ids: List[str],
645    scores: np.ndarray,
646    *,
647    phases: Optional[np.ndarray] = None,
648    phase_tiers: Optional[np.ndarray] = None,
649) -> Dict[str, object]:
650    """Roll per-(scene, phase) scores up to per-scene splits, preserving
651    per-phase tier detail for downstream analysis.
652
653    Each scene's representative score is the **mean of its phase scores**;
654    scenes are sorted by that score and cut into ⌊N/3⌋ Easy + ⌊N/3⌋ Medium
655    + (remainder) Hard. With N=100 the cut is the canonical 33/33/34;
656    with N=99 it becomes 33/33/33.
657
658    Output shape::
659
660        {
661          "splits": {
662            "easy":   [scene_id, …],   # alphabetically sorted in-tier
663            "medium": [...],
664            "hard":   [...]
665          },
666          "scene_detail": {
667            "scene_id": {
668              "scene_tier":  "easy" | "medium" | "hard",
669              "scene_score": <mean of phase scores>,
670              "phase_tiers": ["easy", "medium", "hard"]   # length-3 if phase
671                                                          # info is supplied,
672                                                          # in capture-phase
673                                                          # order (0,1,2)
674              "phase_scores": [0.30, 0.45, 0.28]          # same indexing
675            },
676            ...
677          }
678        }
679
680    The split lists are the consumer-facing deliverable; ``scene_detail``
681    answers "what did each phase look like for this scene?" inline,
682    without requiring readers to cross-reference a separate per-phase
683    difficulty table. Per-phase info is included only when ``phases``
684    (and optionally ``phase_tiers``) are passed.
685    """
686    if len(scene_ids) != scores.size:
687        raise ConfigError(
688            f"length mismatch: {len(scene_ids)} scene_ids vs {scores.size} scores",
689            hint="pass arrays aligned 1:1 across (scene, phase) entries.",
690        )
691    if phases is not None and len(phases) != scores.size:
692        raise ConfigError(
693            f"length mismatch: {len(phases)} phases vs {scores.size} scores",
694            hint="phases must align 1:1 with scene_ids and scores.",
695        )
696    if phase_tiers is not None and len(phase_tiers) != scores.size:
697        raise ConfigError(
698            f"length mismatch: {len(phase_tiers)} phase_tiers vs {scores.size} scores",
699            hint="phase_tiers must align 1:1 with scene_ids and scores.",
700        )
701
702    from collections import defaultdict
703
704    # Aggregate per-scene: mean score plus per-phase detail.
705    per_scene: Dict[str, Dict[int, Tuple[float, Optional[str]]]] = defaultdict(dict)
706    for i, sid in enumerate(scene_ids):
707        ph = int(phases[i]) if phases is not None else i
708        tier_here = str(phase_tiers[i]) if phase_tiers is not None else None
709        per_scene[sid][ph] = (float(scores[i]), tier_here)
710
711    per_scene_mean = {
712        sid: float(np.mean([v[0] for v in info.values()])) for sid, info in per_scene.items()
713    }
714
715    # Sort scenes ascending by mean score; tertile-cut.
716    ordered = sorted(per_scene_mean, key=per_scene_mean.get)
717    n = len(ordered)
718    n_easy = n // 3
719    n_med = n // 3
720    easy = sorted(ordered[:n_easy])
721    medium = sorted(ordered[n_easy : n_easy + n_med])
722    hard = sorted(ordered[n_easy + n_med :])
723    splits = {"easy": easy, "medium": medium, "hard": hard}
724
725    # Reverse map for "which tier does this scene belong to".
726    scene_to_tier = {
727        **{s: "easy" for s in easy},
728        **{s: "medium" for s in medium},
729        **{s: "hard" for s in hard},
730    }
731
732    scene_detail: Dict[str, Dict[str, object]] = {}
733    for sid, info in per_scene.items():
734        # Sort phases by index so output is deterministic (capture order).
735        ordered_phases = sorted(info)
736        phase_scores = [info[p][0] for p in ordered_phases]
737        phase_tiers_out = [info[p][1] for p in ordered_phases if info[p][1] is not None]
738        detail: Dict[str, object] = {
739            "scene_tier": scene_to_tier[sid],
740            "scene_score": per_scene_mean[sid],
741            "phase_scores": phase_scores,
742        }
743        if phase_tiers_out:
744            detail["phase_tiers"] = phase_tiers_out
745        if phases is not None:
746            detail["phase_indices"] = ordered_phases
747        scene_detail[sid] = detail
748
749    return {"splits": splits, "scene_detail": scene_detail}
750
751
752def compute_confidence(
753    pn: np.ndarray,
754    *,
755    n_perturb: int = 1000,
756    seed: int = 0,
757    high_max: float = DEFAULT_CONFIDENCE_HIGH_MAX,
758    medium_max: float = DEFAULT_CONFIDENCE_MEDIUM_MAX,
759) -> Tuple[np.ndarray, np.ndarray]:
760    """Per-row confidence in the uniform mean_pn tier under weight perturbation.
761
762    For each Dirichlet$(\\mathbf{1})$-sampled weight vector, recompute the
763    tertile and check whether each row's tier changed. A row whose tier flips
764    in <``high_max`` of perturbations is "high" confidence; <``medium_max``
765    is "medium"; otherwise "low".
766
767    Returns ``(flip_rate, label)`` arrays, both of length N.
768    """
769    rng = np.random.default_rng(seed)
770    F = pn.shape[1]
771    base_score = pn.mean(axis=1)
772    base_tier = tertile_cut(base_score)
773
774    flips = np.zeros(pn.shape[0], dtype=np.int32)
775    for _ in range(n_perturb):
776        w = rng.dirichlet(np.ones(F))
777        flips += (tertile_cut(pn @ w) != base_tier).astype(np.int32)
778    flip_rate = flips.astype(np.float64) / float(n_perturb)
779
780    label = np.empty(pn.shape[0], dtype=object)
781    label[flip_rate < high_max] = "high"
782    label[(flip_rate >= high_max) & (flip_rate < medium_max)] = "medium"
783    label[flip_rate >= medium_max] = "low"
784    return flip_rate, label