rpx_benchmark.paper_depth_analysis

  1"""Paper-level statistical analysis for the six-metric RPX D1-F cell log."""
  2
  3from __future__ import annotations
  4
  5import csv
  6import hashlib
  7import json
  8import math
  9from pathlib import Path
 10from typing import Any, Mapping, Sequence
 11
 12import numpy as np
 13
 14from .cell_log import merge_cells, write_cells
 15from .exceptions import ConfigError, DatasetError
 16from .metrics.depth_paper import D1_CALIBRATION, PAPER_METRIC_KEYS
 17
 18PHASES = ("clutter", "interaction", "clean")
 19DIFFICULTIES = ("easy", "medium", "hard")
 20METRIC_DIRECTIONS = {
 21    "absrel": "lower",
 22    "rmse": "lower",
 23    "silog": "lower",
 24    "delta1": "higher",
 25    "irmse": "lower",
 26    "fscore_5cm": "higher",
 27}
 28_PHASE_ALIASES = {
 29    "0": "clutter",
 30    "clu": "clutter",
 31    "clutter": "clutter",
 32    "1": "interaction",
 33    "int": "interaction",
 34    "interaction": "interaction",
 35    "2": "clean",
 36    "cln": "clean",
 37    "clean": "clean",
 38}
 39
 40
 41def _phase_contrasts() -> np.ndarray:
 42    return np.array(
 43        [
 44            [-1.0 / math.sqrt(2.0), -1.0 / math.sqrt(6.0)],
 45            [1.0 / math.sqrt(2.0), -1.0 / math.sqrt(6.0)],
 46            [0.0, 2.0 / math.sqrt(6.0)],
 47        ]
 48    )
 49
 50
 51def _residual_ssp(y: np.ndarray, design: np.ndarray) -> tuple[np.ndarray, int]:
 52    rank = int(np.linalg.matrix_rank(design))
 53    df = len(y) - rank
 54    if df <= 0:
 55        raise ConfigError("MANOVA design has no residual degrees of freedom.")
 56    coefficients, *_ = np.linalg.lstsq(design, y, rcond=None)
 57    residual = y - design @ coefficients
 58    return residual.T @ residual, df
 59
 60
 61def _effect_summary(h: np.ndarray, e: np.ndarray, *, df_h: int, df_e: int) -> dict[str, float]:
 62    from scipy.stats import f as f_distribution
 63
 64    h = (h + h.T) / 2.0
 65    e = (e + e.T) / 2.0
 66    response_dim = e.shape[0]
 67    effect_s = min(response_dim, df_h)
 68
 69    # An exactly phase-invariant cube has H=0 and E=0 after subject effects.
 70    # This is a well-defined null effect (Phi=1), even though E has no inverse.
 71    tolerance = np.finfo(np.float64).eps * max(1.0, float(np.linalg.norm(e))) * 1_000
 72    if float(np.linalg.norm(h)) <= tolerance:
 73        p = float(response_dim)
 74        q = float(df_h)
 75        numerator = p * p * q * q - 4.0
 76        denominator = p * p + q * q - 5.0
 77        rao_t = math.sqrt(numerator / denominator) if numerator > 0 and denominator > 0 else 1.0
 78        df1 = response_dim * df_h
 79        df2 = rao_t * (df_e - (response_dim - df_h + 1.0) / 2.0) - (df1 - 2.0) / 2.0
 80        if df2 <= 0:
 81            raise ConfigError("Rao F approximation produced non-positive denominator df.")
 82        return {
 83            "wilks_lambda": 1.0,
 84            "rao_f": 0.0,
 85            "df1": float(df1),
 86            "df2": float(df2),
 87            "p_value": 1.0,
 88            "partial_eta_squared": 0.0,
 89            "phi_wilks": 1.0,
 90            "pillai_trace": 0.0,
 91            "phi_pillai": 1.0,
 92            "pillai_minus_wilks_phi": 0.0,
 93        }
 94    try:
 95        chol = np.linalg.cholesky(e)
 96    except np.linalg.LinAlgError as exc:
 97        raise ConfigError(
 98            "MANOVA residual SSP matrix is singular.",
 99            hint="Check for constant/redundant metrics or incomplete scene-phase cells.",
100        ) from exc
101    transformed = np.linalg.solve(chol, h)
102    transformed = np.linalg.solve(chol, transformed.T).T
103    eigenvalues = np.linalg.eigvalsh((transformed + transformed.T) / 2.0)
104    eigenvalues = np.clip(eigenvalues, 0.0, None)
105    wilks = float(np.prod(1.0 / (1.0 + eigenvalues)))
106    pillai = float(np.sum(eigenvalues / (1.0 + eigenvalues)))
107    phi_wilks = float(wilks ** (1.0 / effect_s))
108    eta_wilks = float(1.0 - phi_wilks)
109    phi_pillai = float(1.0 - pillai / effect_s)
110
111    # Rao's F approximation for Wilks' lambda.
112    p = float(response_dim)
113    q = float(df_h)
114    numerator = p * p * q * q - 4.0
115    denominator = p * p + q * q - 5.0
116    rao_t = math.sqrt(numerator / denominator) if numerator > 0 and denominator > 0 else 1.0
117    df1 = int(response_dim * df_h)
118    df2 = rao_t * (df_e - (response_dim - df_h + 1.0) / 2.0) - (df1 - 2.0) / 2.0
119    if df2 <= 0:
120        raise ConfigError("Rao F approximation produced non-positive denominator df.")
121    lambda_root = wilks ** (1.0 / rao_t)
122    f_value = 0.0 if lambda_root >= 1.0 else ((1.0 - lambda_root) / lambda_root) * (df2 / df1)
123    p_value = float(f_distribution.sf(f_value, df1, df2))
124    return {
125        "wilks_lambda": wilks,
126        "rao_f": float(f_value),
127        "df1": float(df1),
128        "df2": float(df2),
129        "p_value": p_value,
130        "partial_eta_squared": eta_wilks,
131        "phi_wilks": phi_wilks,
132        "pillai_trace": pillai,
133        "phi_pillai": phi_pillai,
134        "pillai_minus_wilks_phi": float(phi_pillai - phi_wilks),
135    }
136
137
138def _nested_effect(
139    y: np.ndarray, reduced: np.ndarray, full: np.ndarray, *, df_h: int
140) -> dict[str, float]:
141    e_reduced, _ = _residual_ssp(y, reduced)
142    e_full, df_e = _residual_ssp(y, full)
143    h = e_reduced - e_full
144    return _effect_summary(h, e_full, df_h=df_h, df_e=df_e)
145
146
147def _repeated_phase_test(cube: np.ndarray) -> dict[str, float]:
148    n_scenes, n_phases, n_metrics = cube.shape
149    if n_phases != 3 or n_metrics != len(PAPER_METRIC_KEYS):
150        raise ConfigError(f"Unexpected MANOVA cube shape: {cube.shape}")
151    y = cube.reshape(n_scenes * n_phases, n_metrics)
152    subjects = np.repeat(np.eye(n_scenes), n_phases, axis=0)
153    phase = np.tile(_phase_contrasts(), (n_scenes, 1))
154    return _nested_effect(y, subjects, np.column_stack((subjects, phase)), df_h=2)
155
156
157def _hotelling_transition(a: np.ndarray, b: np.ndarray) -> dict[str, float]:
158    from scipy.stats import f as f_distribution
159
160    diff = b - a
161    n, k = diff.shape
162    if n <= k:
163        raise ConfigError(f"Hotelling T2 requires scenes > metrics; got {n} <= {k}.")
164    mean = np.mean(diff, axis=0)
165    if np.allclose(diff, 0.0, rtol=0.0, atol=1e-12):
166        return {
167            "hotelling_t_squared": 0.0,
168            "f_value": 0.0,
169            "df1": float(k),
170            "df2": float(n - k),
171            "p_value_raw": 1.0,
172            "wilks_lambda": 1.0,
173            "phi": 1.0,
174        }
175    covariance = np.cov(diff, rowvar=False, ddof=1)
176    try:
177        solved = np.linalg.solve(covariance, mean)
178    except np.linalg.LinAlgError as exc:
179        raise ConfigError("Hotelling covariance matrix is singular.") from exc
180    t_squared = float(n * mean @ solved)
181    f_value = float((n - k) * t_squared / (k * (n - 1)))
182    p_value = float(f_distribution.sf(f_value, k, n - k))
183    wilks = float(1.0 / (1.0 + t_squared / (n - 1)))
184    return {
185        "hotelling_t_squared": t_squared,
186        "f_value": f_value,
187        "df1": float(k),
188        "df2": float(n - k),
189        "p_value_raw": p_value,
190        "wilks_lambda": wilks,
191        "phi": wilks,
192    }
193
194
195def _holm_adjust(values: Sequence[float]) -> list[float]:
196    order = np.argsort(values)
197    adjusted = np.empty(len(values), dtype=np.float64)
198    running = 0.0
199    for rank, index in enumerate(order):
200        candidate = min(1.0, (len(values) - rank) * float(values[index]))
201        running = max(running, candidate)
202        adjusted[index] = running
203    return adjusted.tolist()
204
205
206def _interaction_test(cube: np.ndarray, difficulties: Sequence[str]) -> dict[str, float]:
207    n_scenes, n_phases, n_metrics = cube.shape
208    y = cube.reshape(n_scenes * n_phases, n_metrics)
209    subjects = np.repeat(np.eye(n_scenes), n_phases, axis=0)
210    phase = np.tile(_phase_contrasts(), (n_scenes, 1))
211    difficulty_index = np.array([DIFFICULTIES.index(value) for value in difficulties])
212    difficulty = _phase_contrasts()[difficulty_index]
213    difficulty_long = np.repeat(difficulty, n_phases, axis=0)
214    interaction = (phase[:, :, None] * difficulty_long[:, None, :]).reshape(len(y), 4)
215    reduced = np.column_stack((subjects, phase))
216    full = np.column_stack((reduced, interaction))
217    return _nested_effect(y, reduced, full, df_h=4)
218
219
220def _load_jedi_bounds(path: Path | None) -> tuple[dict[str, Any] | None, str | None]:
221    if path is None:
222        return None, None
223    raw = path.read_bytes()
224    payload = json.loads(raw)
225    if not payload.get("version") or not payload.get("provenance"):
226        raise ConfigError("JEDI bounds require non-empty version and provenance.")
227    epsilon = payload.get("epsilon")
228    if not isinstance(epsilon, (int, float)) or not 0.0 < float(epsilon) <= 1.0:
229        raise ConfigError("JEDI epsilon must be numeric in (0, 1].")
230    metrics = payload.get("metrics")
231    if not isinstance(metrics, dict) or set(metrics) != set(PAPER_METRIC_KEYS):
232        raise ConfigError(f"JEDI bounds must define exactly: {list(PAPER_METRIC_KEYS)}")
233    for key in PAPER_METRIC_KEYS:
234        spec = metrics[key]
235        if spec.get("direction") != METRIC_DIRECTIONS[key]:
236            raise ConfigError(f"JEDI direction mismatch for {key}.")
237        if not all(isinstance(spec.get(field), (int, float)) for field in ("best", "worst")):
238            raise ConfigError(f"JEDI {key} requires numeric best and worst bounds.")
239        best = float(spec["best"])
240        worst = float(spec["worst"])
241        if not np.isfinite([best, worst]).all():
242            raise ConfigError(f"JEDI {key} bounds must be finite.")
243        if (spec["direction"] == "lower" and best >= worst) or (
244            spec["direction"] == "higher" and best <= worst
245        ):
246            raise ConfigError(f"JEDI {key} best/worst bounds contradict its direction.")
247    return payload, hashlib.sha256(raw).hexdigest()
248
249
250def _jedi_score(values: Mapping[str, float], bounds: Mapping[str, Any]) -> dict[str, Any]:
251    epsilon = float(bounds["epsilon"])
252    desirabilities: dict[str, float] = {}
253    for key in PAPER_METRIC_KEYS:
254        spec = bounds["metrics"][key]
255        value = float(values[key])
256        desirability = (value - float(spec["worst"])) / (float(spec["best"]) - float(spec["worst"]))
257        desirabilities[key] = float(np.clip(desirability, 0.0, 1.0))
258    score = float(
259        np.exp(np.mean(np.log([max(value, epsilon) for value in desirabilities.values()])))
260    )
261    return {"score": score, "desirabilities": desirabilities}
262
263
264def analyze_depth_cells(
265    paths: Sequence[str | Path],
266    *,
267    jedi_bounds_path: str | Path | None = None,
268    expected_scenes: int = 100,
269    provenance: Mapping[str, Any] | None = None,
270) -> tuple[dict[str, Any], list[dict[str, Any]]]:
271    """Validate, merge and analyze three RPX D1-F split cell logs."""
272    rows = merge_cells(Path(path) for path in paths)
273    rows = [
274        row
275        for row in rows
276        if row.get("task") == "monocular_depth" and int(row.get("frame_budget") or 0) == 0
277    ]
278    if not rows:
279        raise DatasetError("No monocular_depth cells found in the supplied logs.")
280    model_names = {str(row.get("model_name")) for row in rows}
281    if len(model_names) != 1 or model_names == {"None"} or model_names == {""}:
282        raise DatasetError(f"Paper analysis requires exactly one model; got {sorted(model_names)}")
283
284    by_key: dict[tuple[str, str], dict[str, Any]] = {}
285    for row in rows:
286        scene = str(row.get("scene_id") or "")
287        phase_raw = str(row.get("phase") or "").lower()
288        phase = _PHASE_ALIASES.get(phase_raw)
289        if not scene or phase is None:
290            raise DatasetError(f"Invalid scene/phase cell: {scene!r}/{phase_raw!r}")
291        key = (scene, phase)
292        if key in by_key:
293            raise DatasetError(f"Duplicate scene-phase cell: {key}")
294        normalized = dict(row)
295        normalized["phase"] = phase
296        by_key[key] = normalized
297
298    scenes = sorted({scene for scene, _phase in by_key})
299    if len(scenes) != expected_scenes:
300        raise DatasetError(f"Expected {expected_scenes} scenes; found {len(scenes)}.")
301    if len(by_key) != expected_scenes * 3:
302        raise DatasetError(f"Expected {expected_scenes * 3} cells; found {len(by_key)}.")
303
304    raw_cube = np.empty((len(scenes), 3, len(PAPER_METRIC_KEYS)), dtype=np.float64)
305    difficulties: list[str] = []
306    normalized_rows: list[dict[str, Any]] = []
307    for scene_index, scene in enumerate(scenes):
308        scene_difficulties = set()
309        for phase_index, phase in enumerate(PHASES):
310            if (scene, phase) not in by_key:
311                raise DatasetError(f"Missing scene-phase cell: {scene}/{phase}")
312            row = by_key[(scene, phase)]
313            difficulty = str(row.get("difficulty") or "").lower()
314            if difficulty not in DIFFICULTIES:
315                raise DatasetError(f"Invalid difficulty for {scene}/{phase}: {difficulty!r}")
316            scene_difficulties.add(difficulty)
317            for metric_index, metric in enumerate(PAPER_METRIC_KEYS):
318                value = row.get(f"metric:{metric}")
319                if not isinstance(value, (int, float)) or not np.isfinite(value):
320                    raise DatasetError(f"Missing/non-finite {metric} for {scene}/{phase}.")
321                raw_cube[scene_index, phase_index, metric_index] = float(value)
322            normalized_rows.append(row)
323        if len(scene_difficulties) != 1:
324            raise DatasetError(f"Difficulty changes across phases for scene {scene}.")
325        difficulties.append(scene_difficulties.pop())
326
327    oriented = raw_cube.copy()
328    for index, metric in enumerate(PAPER_METRIC_KEYS):
329        if METRIC_DIRECTIONS[metric] == "lower":
330            oriented[:, :, index] *= -1.0
331    flat = oriented.reshape(-1, len(PAPER_METRIC_KEYS))
332    std = np.std(flat, axis=0, ddof=0)
333    if np.any(std <= 0) or not np.isfinite(std).all():
334        raise DatasetError("Paper analysis requires six non-constant finite metrics.")
335    standardized = ((flat - np.mean(flat, axis=0)) / std).reshape(oriented.shape)
336
337    overall = _repeated_phase_test(standardized)
338    transitions: list[dict[str, Any]] = []
339    for label, first, second in (
340        ("clutter_to_interaction", 0, 1),
341        ("interaction_to_clean", 1, 2),
342        ("clutter_to_clean", 0, 2),
343    ):
344        transitions.append(
345            {
346                "transition": label,
347                **_hotelling_transition(standardized[:, first, :], standardized[:, second, :]),
348            }
349        )
350    adjusted = _holm_adjust([item["p_value_raw"] for item in transitions])
351    for item, value in zip(transitions, adjusted, strict=True):
352        item["p_value_holm"] = value
353
354    per_tier: dict[str, Any] = {}
355    difficulty_array = np.asarray(difficulties)
356    for tier in DIFFICULTIES:
357        subset = standardized[difficulty_array == tier]
358        per_tier[tier] = {"n_scenes": int(len(subset)), **_repeated_phase_test(subset)}
359
360    phase_means = {
361        phase: {
362            metric: float(np.mean(raw_cube[:, phase_index, metric_index]))
363            for metric_index, metric in enumerate(PAPER_METRIC_KEYS)
364        }
365        for phase_index, phase in enumerate(PHASES)
366    }
367    overall_means = {
368        metric: float(np.mean(raw_cube[:, :, metric_index]))
369        for metric_index, metric in enumerate(PAPER_METRIC_KEYS)
370    }
371
372    bounds_path = Path(jedi_bounds_path) if jedi_bounds_path else None
373    bounds, bounds_sha = _load_jedi_bounds(bounds_path)
374    if bounds is None:
375        jedi: dict[str, Any] = {"status": "bounds_missing"}
376    else:
377        jedi = {
378            "status": "computed",
379            "bounds_version": bounds["version"],
380            "bounds_sha256": bounds_sha,
381            "overall": _jedi_score(overall_means, bounds),
382            "per_phase": {
383                phase: _jedi_score(values, bounds) for phase, values in phase_means.items()
384            },
385        }
386
387    analysis = {
388        "schema_version": "rpx-d1f-paper-analysis-v1",
389        "model_name": next(iter(model_names)),
390        "metrics": list(PAPER_METRIC_KEYS),
391        "metric_directions": METRIC_DIRECTIONS,
392        "n_scenes": len(scenes),
393        "n_cells": len(by_key),
394        "calibration": D1_CALIBRATION.to_dict(),
395        "overall_phase_manova": overall,
396        "transitions": transitions,
397        "per_tier": per_tier,
398        "phase_by_difficulty": _interaction_test(standardized, difficulties),
399        "metric_means": {"overall": overall_means, "per_phase": phase_means},
400        "jedi": jedi,
401        "jedi_status": jedi["status"],
402        "provenance": dict(provenance or {}),
403    }
404    return analysis, normalized_rows
405
406
407def write_depth_analysis(
408    analysis: Mapping[str, Any],
409    cells: Sequence[Mapping[str, Any]],
410    output_dir: str | Path,
411) -> dict[str, Path]:
412    out = Path(output_dir)
413    out.mkdir(parents=True, exist_ok=True)
414    combined = out / "combined_cells.parquet"
415    write_cells(cells, combined, format="parquet")
416    json_path = out / "paper_analysis.json"
417    json_path.write_text(json.dumps(analysis, indent=2, sort_keys=True) + "\n")
418
419    table_path = out / "paper_table.csv"
420    flat: dict[str, Any] = {
421        "model_name": analysis["model_name"],
422        "phi": analysis["overall_phase_manova"]["phi_wilks"],
423        "wilks_lambda": analysis["overall_phase_manova"]["wilks_lambda"],
424        "rao_f": analysis["overall_phase_manova"]["rao_f"],
425        "p_value": analysis["overall_phase_manova"]["p_value"],
426        "partial_eta_squared": analysis["overall_phase_manova"]["partial_eta_squared"],
427        "phi_pillai": analysis["overall_phase_manova"]["phi_pillai"],
428        "wilks_pillai_gap": analysis["overall_phase_manova"]["pillai_minus_wilks_phi"],
429        "jedi_status": analysis["jedi"]["status"],
430    }
431    if analysis["jedi"]["status"] == "computed":
432        flat["jedi"] = analysis["jedi"]["overall"]["score"]
433    for phase, values in analysis["metric_means"]["per_phase"].items():
434        for metric, value in values.items():
435            flat[f"{metric}_{phase}"] = value
436    for item in analysis["transitions"]:
437        flat[f"phi_{item['transition']}"] = item["phi"]
438        flat[f"hotelling_t2_{item['transition']}"] = item["hotelling_t_squared"]
439        flat[f"p_raw_{item['transition']}"] = item["p_value_raw"]
440        flat[f"p_holm_{item['transition']}"] = item["p_value_holm"]
441    for tier, values in analysis["per_tier"].items():
442        for key in ("n_scenes", "wilks_lambda", "phi_wilks", "p_value", "partial_eta_squared"):
443            flat[f"{key}_{tier}"] = values[key]
444    for key in ("wilks_lambda", "rao_f", "p_value", "partial_eta_squared"):
445        flat[f"phase_by_difficulty_{key}"] = analysis["phase_by_difficulty"][key]
446    with table_path.open("w", newline="") as handle:
447        writer = csv.DictWriter(handle, fieldnames=list(flat))
448        writer.writeheader()
449        writer.writerow(flat)
450
451    markdown_path = out / "paper_analysis.md"
452    overall = analysis["overall_phase_manova"]
453    lines = [
454        "# RPX D1-F Paper Analysis",
455        "",
456        f"- Model: `{analysis['model_name']}`",
457        f"- Complete cells: {analysis['n_cells']} ({analysis['n_scenes']} scenes x 3 phases)",
458        f"- Wilks Lambda: {overall['wilks_lambda']:.6g}",
459        f"- Phi: {overall['phi_wilks']:.6g}",
460        f"- Rao F({overall['df1']:.0f}, {overall['df2']:.3f}): "
461        f"{overall['rao_f']:.6g}; p={overall['p_value']:.6g}",
462        f"- Partial eta squared: {overall['partial_eta_squared']:.6g}",
463        f"- Pillai Phi: {overall['phi_pillai']:.6g}; gap={overall['pillai_minus_wilks_phi']:.6g}",
464        f"- JEDI: `{analysis['jedi']['status']}`",
465        "",
466        "## Phase transitions",
467        "",
468        "| Transition | Hotelling T2 | Wilks Phi | Raw p | Holm p |",
469        "|---|---:|---:|---:|---:|",
470    ]
471    for item in analysis["transitions"]:
472        lines.append(
473            f"| {item['transition']} | {item['hotelling_t_squared']:.6g} | "
474            f"{item['phi']:.6g} | {item['p_value_raw']:.6g} | "
475            f"{item['p_value_holm']:.6g} |"
476        )
477    lines.extend(
478        [
479            "",
480            "## Difficulty tiers",
481            "",
482            "| Tier | Scenes | Wilks Lambda | Phi | p | Partial eta squared |",
483            "|---|---:|---:|---:|---:|---:|",
484        ]
485    )
486    for tier in DIFFICULTIES:
487        values = analysis["per_tier"][tier]
488        lines.append(
489            f"| {tier} | {values['n_scenes']} | {values['wilks_lambda']:.6g} | "
490            f"{values['phi_wilks']:.6g} | {values['p_value']:.6g} | "
491            f"{values['partial_eta_squared']:.6g} |"
492        )
493    interaction = analysis["phase_by_difficulty"]
494    lines.extend(
495        [
496            "",
497            "## Phase by difficulty",
498            "",
499            f"Wilks Lambda={interaction['wilks_lambda']:.6g}; "
500            f"Rao F({interaction['df1']:.0f}, {interaction['df2']:.3f})="
501            f"{interaction['rao_f']:.6g}; p={interaction['p_value']:.6g}; "
502            f"partial eta squared={interaction['partial_eta_squared']:.6g}.",
503            "",
504            "## Phase means",
505            "",
506            "| Phase | " + " | ".join(PAPER_METRIC_KEYS) + " |",
507            "|---|" + "---:|" * len(PAPER_METRIC_KEYS),
508        ]
509    )
510    for phase in PHASES:
511        values = analysis["metric_means"]["per_phase"][phase]
512        lines.append(
513            "| "
514            + phase
515            + " | "
516            + " | ".join(f"{values[key]:.6g}" for key in PAPER_METRIC_KEYS)
517            + " |"
518        )
519    markdown_path.write_text("\n".join(lines) + "\n")
520    return {
521        "combined_cells": combined,
522        "analysis_json": json_path,
523        "analysis_markdown": markdown_path,
524        "paper_table": table_path,
525    }