operator_graph_validation.py

#!/usr/bin/env python3
"""
Operator-graph validation for NeRArch-Sim.

Validates a captured/transformed neural-rendering operator graph along three axes:

  1. Structural  - DAG acyclicity, taxonomy-stage coverage/order, renderer source
                   wiring, and a non-degenerate op-type histogram. Works with no
                   capture data (operator graph only).
  2. Functional  - Replays the captured per-sample field outputs into an image and
                   compares (PSNR) against the nerfstudio-rendered reference panel.
                   Two modes:
                     * checkpoint     : trusts captured weights + per-sample colors.
                     * force_recompute: ignores captured weights and recomputes them
                                        from captured density + sample geometry
                                        (deltas), proving the sampling->field->
                                        blending chain (anti-cheat).
  3. Ablation    - Zeroes/perturbs each stage of the forced-recompute pipeline and
                   requires PSNR to drop by >= margin, proving every stage materially
                   contributes (a bypassed stage => FAIL).

Verdict: PASS iff structural PASS AND force-recompute PSNR >= psnr_threshold AND
every ablation drop >= ablation_margin.

Scope: NeRF (volumetric) replay. 3DGS/splatfacto alpha-blend replay is a documented
follow-up; for such captures functional validation reports "unsupported" and the
verdict falls back to structural-only.
"""

from __future__ import annotations

import json
import math
from pathlib import Path
from typing import Any, Dict, List, Optional, Tuple

PSNR_THRESHOLD_DEFAULT = 35.0
ABLATION_MARGIN_DEFAULT = 10.0
GT_PARITY_MARGIN_DEFAULT = 1.0

# Taxonomy stage order used for structural ordering checks.
_STAGE_ORDER = ["SAMPLING", "ENCODING", "FIELD_COMPUTATION", "BLENDING"]
_STAGE_RANK = {s: i for i, s in enumerate(_STAGE_ORDER)}

# Stages whose ablation we exercise in the forced-recompute pipeline.
_ABLATIONS = ["sampling", "field_density", "field_color", "weights"]


# --------------------------------------------------------------------------- #
# Small helpers
# --------------------------------------------------------------------------- #
def _psnr(a, b) -> float:
    import numpy as np

    a = a.astype(np.float32)
    b = b.astype(np.float32)
    mse = float(np.mean((a - b) ** 2))
    if mse <= 1e-12:
        return 99.0
    return float(20.0 * math.log10(255.0) - 10.0 * math.log10(mse))


def _taxonomy_stage(op_obj) -> str:
    """Map an /Operators node to a unified taxonomy stage by class name."""
    cn = type(op_obj).__name__
    if any(t in cn for t in ["Sampler", "Sample", "FrustumCulling", "FrustrumCulling", "Culling"]):
        return "SAMPLING"
    if any(t in cn for t in ["Encoding", "Encoder", "Hash", "RFF", "Positional", "Fourier", "SphericalHarm"]):
        return "ENCODING"
    if any(t in cn for t in ["Render", "Blend", "Volume", "RGB", "Alpha", "Composite"]):
        return "BLENDING"
    if any(t in cn for t in ["MLP", "Network", "Field", "Computation", "Density", "Color"]):
        return "FIELD_COMPUTATION"
    return "FIELD_COMPUTATION"


# --------------------------------------------------------------------------- #
# Structural validation
# --------------------------------------------------------------------------- #
def validate_structure(operators_graph) -> Dict[str, Any]:
    """Deterministic structural checks over the transformed operator graph."""
    nodes = list(operators_graph.nodes)
    children_map: Dict[Any, List[Any]] = {n: list(getattr(n, "children", []) or []) for n in nodes}

    indeg: Dict[Any, int] = {n: 0 for n in nodes}
    parents: Dict[Any, List[Any]] = {n: [] for n in nodes}
    for n in nodes:
        for c in children_map[n]:
            if c in indeg:
                indeg[c] += 1
                parents[c].append(n)

    # --- DAG acyclicity (Kahn topological sort) ---
    from collections import deque

    indeg_work = dict(indeg)
    q = deque([n for n in nodes if indeg_work[n] == 0])
    topo: List[Any] = []
    while q:
        n = q.popleft()
        topo.append(n)
        for c in children_map[n]:
            if c in indeg_work:
                indeg_work[c] -= 1
                if indeg_work[c] == 0:
                    q.append(c)
    is_dag = len(topo) == len(nodes)

    # --- Taxonomy histogram + stage coverage ---
    stage_hist: Dict[str, int] = {}
    type_hist: Dict[str, int] = {}
    for n in nodes:
        s = _taxonomy_stage(n)
        stage_hist[s] = stage_hist.get(s, 0) + 1
        cn = type(n).__name__
        type_hist[cn] = type_hist.get(cn, 0) + 1

    present_stages = {s for s, c in stage_hist.items() if c > 0}

    # --- Stage-order edges (forward vs backward) ---
    forward_edges = 0
    backward_edges = 0
    for u in nodes:
        ru = _STAGE_RANK.get(_taxonomy_stage(u), 99)
        for v in children_map[u]:
            rv = _STAGE_RANK.get(_taxonomy_stage(v), 99)
            if rv >= ru:
                forward_edges += 1
            else:
                backward_edges += 1
    total_edges = forward_edges + backward_edges
    # Backward edges are legitimate for coarse->fine resampling, but should be a
    # minority of all edges.
    forward_ratio = (forward_edges / total_edges) if total_edges else 1.0

    # --- Renderer source wiring: every BLENDING node must have a predecessor ---
    blending_nodes = [n for n in nodes if _taxonomy_stage(n) == "BLENDING"]
    blending_with_source = [n for n in blending_nodes if len(parents.get(n, [])) > 0]
    sinks = [n for n in nodes if len(children_map[n]) == 0]
    blending_sinks = [n for n in sinks if _taxonomy_stage(n) == "BLENDING"]

    checks: List[Dict[str, Any]] = []

    def _add(name: str, passed: bool, detail: str) -> None:
        checks.append({"name": name, "pass": bool(passed), "detail": detail})

    _add("dag_acyclic", is_dag,
         f"topological order covered {len(topo)}/{len(nodes)} nodes")
    _add("stages_present", present_stages >= {"SAMPLING", "FIELD_COMPUTATION", "BLENDING"},
         f"present stages: {sorted(present_stages)}")
    _add("has_blending_sink", len(blending_sinks) > 0,
         f"{len(blending_sinks)} BLENDING sink(s)")
    _add("blending_sources_wired",
         len(blending_nodes) > 0 and len(blending_with_source) == len(blending_nodes),
         f"{len(blending_with_source)}/{len(blending_nodes)} BLENDING nodes have predecessors")
    _add("stage_order_forward_majority", forward_ratio >= 0.5,
         f"forward edges {forward_edges}/{total_edges} ({forward_ratio:.2f})")
    _add("op_type_histogram_nondegenerate", len(present_stages) >= 2 and len(nodes) >= 4,
         f"{len(nodes)} nodes across {len(present_stages)} stages, {len(type_hist)} op types")

    structural_pass = all(c["pass"] for c in checks)

    return {
        "pass": bool(structural_pass),
        "checks": checks,
        "total_nodes": len(nodes),
        "total_edges": total_edges,
        "stage_histogram": stage_hist,
        "op_type_histogram": type_hist,
        "forward_edges": forward_edges,
        "backward_edges": backward_edges,
        "blending_nodes": len(blending_nodes),
        "blending_sinks": len(blending_sinks),
    }


# --------------------------------------------------------------------------- #
# Capture loading + functional reconstruction
# --------------------------------------------------------------------------- #
def _resolve_chunk_file(raw_path: str, capture_dir: Path) -> Path:
    p = Path(raw_path)
    if p.exists():
        return p
    cand = capture_dir / p.name
    return cand if cand.exists() else p


def load_capture(dag_dir: Path) -> Optional[Dict[str, Any]]:
    """Load the render_capture manifest and group chunks by kind (fine stage)."""
    manifest_path = dag_dir / "render_capture" / "manifest.json"
    if not manifest_path.exists():
        return None
    manifest = json.loads(manifest_path.read_text())
    chunks = manifest.get("chunks", []) or []

    def grab(kind: str) -> List[Dict[str, Any]]:
        items = [c for c in chunks if str(c.get("kind", "")) == kind]
        fine = [c for c in items if str(c.get("stage", "")).lower() == "fine"]
        items = fine if fine else items
        return sorted(items, key=lambda c: int(c.get("index", 0)))

    return {
        "manifest_path": manifest_path,
        "capture_dir": manifest_path.parent,
        "sample": grab("sample_source"),
        "density": grab("density_field_output"),
        "rgb": grab("rgb_field_output"),
        "weights": grab("weights_output"),
        "chunks_total": len(chunks),
    }


def capture_supported(bundle: Optional[Dict[str, Any]]) -> bool:
    """NeRF volumetric replay requires density+rgb+weights+sample chunks."""
    if not bundle:
        return False
    return bool(bundle["density"] and bundle["rgb"] and bundle["weights"] and bundle["sample"])


def _reconstruct_all(bundle: Dict[str, Any], modes: List[str]) -> Dict[str, Any]:
    """
    Single pass over capture chunks producing the per-ray object color term and
    accumulated opacity for every requested reconstruction mode.

    Returns dict with per-mode {"obj": (R,3), "acc": (R,)} plus "n_rays".
    Composition with a background colour is applied later by the caller.
    """
    import numpy as np

    cap = bundle["capture_dir"]

    def load(c, key="value"):
        f = _resolve_chunk_file(c["file"], cap)
        return np.load(f)[key].astype(np.float32)

    def load_sample(c):
        f = _resolve_chunk_file(c["file"], cap)
        d = np.load(f)
        return d["starts"].astype(np.float32), d["ends"].astype(np.float32)

    acc_obj: Dict[str, List] = {m: [] for m in modes}
    acc_op: Dict[str, List] = {m: [] for m in modes}

    groups = list(zip(bundle["sample"], bundle["density"], bundle["rgb"], bundle["weights"]))
    for sc, dc, rc, wc in groups:
        starts, ends = load_sample(sc)
        starts = starts[..., 0]
        ends = ends[..., 0]
        sigma = load(dc)
        if sigma.ndim == 3:
            sigma = sigma[..., 0]
        color = load(rc)[..., :3]
        w_cap = load(wc)
        if w_cap.ndim == 3:
            w_cap = w_cap[..., 0]

        n = min(starts.shape[1], ends.shape[1], sigma.shape[1], color.shape[1], w_cap.shape[1])
        starts, ends = starts[:, :n], ends[:, :n]
        sigma, color, w_cap = sigma[:, :n], color[:, :n], w_cap[:, :n]
        deltas0 = ends - starts

        for mode in modes:
            deltas = deltas0
            sg = sigma
            col = color
            if mode == "abl_sampling":
                deltas = np.full_like(deltas0, float(deltas0.mean()) if deltas0.size else 0.0)
            elif mode == "abl_field_density":
                sg = np.zeros_like(sigma)
            elif mode == "abl_field_color":
                col = np.zeros_like(color)

            if mode == "checkpoint":
                w = w_cap
            else:
                alpha = 1.0 - np.exp(-sg * deltas)
                trans = np.cumprod(1.0 - alpha + 1e-10, axis=1)
                trans = np.concatenate([np.ones((trans.shape[0], 1), np.float32), trans[:, :-1]], axis=1)
                w = alpha * trans
                if mode == "abl_weights":
                    w = np.full_like(w, 1.0 / max(n, 1))

            acc_obj[mode].append(np.sum(w[..., None] * col, axis=1))
            acc_op[mode].append(np.sum(w, axis=1))

    out: Dict[str, Any] = {}
    n_rays = 0
    for mode in modes:
        obj = np.concatenate(acc_obj[mode], axis=0)
        op = np.concatenate(acc_op[mode], axis=0)
        n_rays = obj.shape[0]
        out[mode] = {"obj": obj, "acc": op}
    out["n_rays"] = n_rays
    return out


def _compose(obj, acc, bg: float, h: int, w: int):
    import numpy as np

    need = h * w
    composed = obj + (1.0 - acc)[:, None] * bg
    if composed.shape[0] < need:
        pad = np.ones((need - composed.shape[0], 3), np.float32) * bg
        composed = np.concatenate([composed, pad], axis=0)
    img = np.clip(composed[:need] * 255.0, 0, 255).reshape(h, w, 3).astype(np.uint8)
    return img


def _find_reference(dag_dir: Path) -> Optional[Path]:
    candidates = sorted(dag_dir.rglob("eval_img*.png"))
    if candidates:
        return candidates[0]
    # Fallback: any non-artifact rendered image.
    for p in sorted(dag_dir.rglob("*.png")):
        n = p.name.lower()
        if any(t in n for t in ["operator_graph", "validation", "replay", "ablation", "transformation"]):
            continue
        return p
    return None


def run_functional(
    bundle: Dict[str, Any],
    dag_dir: Path,
    output_prefix: Path,
    psnr_threshold: float,
    ablation_margin: float,
    do_ablation: bool = True,
    gt_parity_margin: float = GT_PARITY_MARGIN_DEFAULT,
    check_gt_parity: bool = True,
) -> Dict[str, Any]:
    """Replay the captured graph, compare to the render + GT panels, run ablations.

    The nerfstudio reference is a montage ``cat([image, rgb_coarse, rgb_fine])``;
    panel 0 is ground truth and the panel best matching our reconstruction is the
    nerfstudio render. We therefore report three PSNRs:
      * replay-vs-render  -> extraction fidelity (gates the verdict),
      * replay-vs-GT      -> our render quality (informational),
      * nerfstudio-vs-GT  -> nerfstudio's own render quality (reference),
    plus a GT-parity check (|replay-vs-GT - nerfstudio-vs-GT| <= margin) that
    confirms our replay does not lose quality relative to nerfstudio.
    """
    import numpy as np
    from PIL import Image

    result: Dict[str, Any] = {
        "supported": True,
        "reference_image": None,
        "reference_panel_index": None,  # alias of render_panel_index (back-compat)
        "render_panel_index": None,
        "gt_panel_index": None,
        "background": None,
        "n_rays": 0,
        "panel_psnr_checkpoint": [],
        "psnr_checkpoint": None,
        "psnr_force_recompute": None,    # alias of psnr_replay_vs_render (back-compat)
        "psnr_replay_vs_render": None,
        "psnr_replay_vs_gt": None,
        "psnr_nerfstudio_vs_gt": None,
        "gt_parity_delta": None,
        "gt_parity_pass": None,
        "gt_parity_checked": False,
        "ablations": [],
        "replay_image": None,
        "reference_panel_image": None,   # alias of nerfstudio_render_image (back-compat)
        "nerfstudio_render_image": None,
        "ground_truth_image": None,
        "error": None,
    }

    ref_path = _find_reference(dag_dir)
    if ref_path is None:
        result["supported"] = False
        result["error"] = "no reference image (eval_img*.png) found"
        return result
    result["reference_image"] = str(ref_path)
    ref = np.array(Image.open(ref_path).convert("RGB")).astype(np.float32)
    ref_h, ref_w = ref.shape[:2]

    modes = ["checkpoint", "recompute"]
    if do_ablation:
        modes += ["abl_" + a for a in _ABLATIONS]
    recon = _reconstruct_all(bundle, modes)
    n_rays = recon["n_rays"]
    result["n_rays"] = int(n_rays)

    # The reference is a horizontal montage of equally sized panels. Derive the
    # panel width from the replayed ray count (one panel = one rendered image of
    # ref_h rows) rather than assuming square panels -- montage panels are only
    # square for synthetic 800x800 scenes; real 360 scenes are e.g. 779x520.
    h = ref_h
    render_w = max(1, n_rays // ref_h)
    panel_w = min(render_w, ref_w)
    n_panels = max(1, int(round(ref_w / panel_w)))
    w = panel_w

    def _panel(idx: int):
        return ref[:, idx * panel_w:(idx + 1) * panel_w, :][:, :w, :]

    # The nerfstudio montage is cat([gt, (coarse,) fine]); ground truth is panel 0.
    # We select the *render* panel as the one our reconstruction best matches,
    # restricted to panels 1..n-1 so it is never confused with ground truth.
    has_gt = n_panels >= 2
    gt_panel_idx = 0 if has_gt else None
    render_candidates = range(1, n_panels) if has_gt else range(n_panels)

    # Select the render panel + background by best checkpoint PSNR.
    ck = recon["checkpoint"]
    best = (-1.0, next(iter(render_candidates)), 1.0)  # (psnr, panel_idx, bg)
    panel_scores: List[Dict[str, Any]] = []
    for bg in (1.0, 0.0):
        img = _compose(ck["obj"], ck["acc"], bg, h, w)
        for p in range(n_panels):
            val = _psnr(img, _panel(p))
            panel_scores.append({"panel": p, "bg": "white" if bg == 1.0 else "black", "psnr": round(val, 3)})
            if p in render_candidates and val > best[0]:
                best = (val, p, bg)
    _, render_idx, bg = best
    result["panel_psnr_checkpoint"] = panel_scores
    result["render_panel_index"] = int(render_idx)
    result["reference_panel_index"] = int(render_idx)
    result["gt_panel_index"] = (int(gt_panel_idx) if gt_panel_idx is not None else None)
    result["background"] = "white" if bg == 1.0 else "black"

    render_panel = _panel(render_idx)
    gt_panel = _panel(gt_panel_idx) if has_gt else None

    ck_img = _compose(ck["obj"], ck["acc"], bg, h, w)
    rc = recon["recompute"]
    rc_img = _compose(rc["obj"], rc["acc"], bg, h, w)

    # Fidelity: our replay vs nerfstudio's render (gates the verdict).
    result["psnr_checkpoint"] = round(_psnr(ck_img, render_panel), 3)
    psnr_replay_vs_render = round(_psnr(rc_img, render_panel), 3)
    result["psnr_replay_vs_render"] = psnr_replay_vs_render
    result["psnr_force_recompute"] = psnr_replay_vs_render
    base_psnr = psnr_replay_vs_render

    # Quality + GT parity (only when a ground-truth panel exists).
    if has_gt:
        psnr_replay_vs_gt = round(_psnr(rc_img, gt_panel), 3)
        psnr_nerf_vs_gt = round(_psnr(render_panel, gt_panel), 3)
        result["psnr_replay_vs_gt"] = psnr_replay_vs_gt
        result["psnr_nerfstudio_vs_gt"] = psnr_nerf_vs_gt
        delta = round(abs(psnr_replay_vs_gt - psnr_nerf_vs_gt), 3)
        result["gt_parity_delta"] = delta
        if check_gt_parity:
            result["gt_parity_checked"] = True
            result["gt_parity_pass"] = bool(delta <= gt_parity_margin)
        else:
            result["gt_parity_pass"] = None

    # Save the 3-way image strip: ground truth | nerfstudio render | replay.
    replay_path = Path(f"{output_prefix}_replay.png")
    Image.fromarray(rc_img).save(replay_path)
    result["replay_image"] = str(replay_path)
    render_path = Path(f"{output_prefix}_nerfstudio_render.png")
    Image.fromarray(render_panel.astype(np.uint8)).save(render_path)
    result["nerfstudio_render_image"] = str(render_path)
    result["reference_panel_image"] = str(render_path)
    if has_gt:
        gt_path = Path(f"{output_prefix}_ground_truth.png")
        Image.fromarray(gt_panel.astype(np.uint8)).save(gt_path)
        result["ground_truth_image"] = str(gt_path)

    if do_ablation:
        for a in _ABLATIONS:
            ab = recon["abl_" + a]
            ab_img = _compose(ab["obj"], ab["acc"], bg, h, w)
            ab_psnr = round(_psnr(ab_img, render_panel), 3)
            drop = round(base_psnr - ab_psnr, 3)
            img_path = Path(f"{output_prefix}_ablation_{a}.png")
            Image.fromarray(ab_img).save(img_path)
            result["ablations"].append({
                "stage": a,
                "psnr": ab_psnr,
                "psnr_drop": drop,
                "pass": bool(drop >= ablation_margin),
                "image": str(img_path),
            })

    return result


# --------------------------------------------------------------------------- #
# Orchestration
# --------------------------------------------------------------------------- #
def validate_operator_graph(
    dag_path: str,
    output_prefix: Optional[str] = None,
    psnr_threshold: float = PSNR_THRESHOLD_DEFAULT,
    ablation_margin: float = ABLATION_MARGIN_DEFAULT,
    do_ablation: bool = True,
    gt_parity_margin: float = GT_PARITY_MARGIN_DEFAULT,
    check_gt_parity: bool = True,
    verbose: bool = True,
) -> Dict[str, Any]:
    """
    Full operator-graph validation entry point.

    Loads + transforms the traced DAG, runs structural checks, and (when capture
    data is available) functional replay + ablation. Writes a JSON + Markdown
    report next to the DAG and returns the report dict (with overall verdict).
    """
    dag_file = Path(dag_path).resolve()
    dag_dir = dag_file.parent
    if output_prefix is None:
        output_prefix = str(dag_dir / "operator_graph_validation_report")
    out_prefix = Path(output_prefix)

    report: Dict[str, Any] = {
        "dag_path": str(dag_file),
        "thresholds": {
            "psnr_threshold_db": psnr_threshold,
            "ablation_margin_db": ablation_margin,
            "gt_parity_margin_db": gt_parity_margin if check_gt_parity else None,
        },
        "structural": None,
        "functional": None,
        "verdict": "FAIL",
        "verdict_reasons": [],
    }

    # --- Transform DAG -> operator graph for structural checks ---
    try:
        try:
            from dag_to_operators_integration import DAGToOperatorsIntegration
        except ImportError:
            from Instrumentation.dag_to_operators_integration import DAGToOperatorsIntegration
        import pickle
        import networkx as nx  # noqa: F401

        with open(dag_file, "rb") as f:
            dag_data = pickle.load(f)
        if hasattr(dag_data, "nodes") and not isinstance(dag_data, dict):
            dict_dag = {"nodes": {}, "edges": list(dag_data.edges())}
            for node_id, node_data in dag_data.nodes(data=True):
                nd = dict(node_data)
                nd["function_name"] = str(node_id)
                dict_dag["nodes"][node_id] = nd
            dag_data = dict_dag
        integration = DAGToOperatorsIntegration()
        operators_graph, _ = integration.transform_dag_to_operators(dag_data)
        report["structural"] = validate_structure(operators_graph)
        if verbose:
            print(f"   [structural] {'PASS' if report['structural']['pass'] else 'FAIL'} "
                  f"({report['structural']['total_nodes']} nodes)")
    except Exception as e:  # pragma: no cover - defensive
        report["structural"] = {"pass": False, "error": str(e), "checks": []}
        if verbose:
            print(f"   [structural] ERROR: {e}")

    # --- Functional replay (requires capture) ---
    bundle = load_capture(dag_dir)
    if not capture_supported(bundle):
        report["functional"] = {
            "supported": False,
            "error": "no NeRF volumetric capture (render_capture/manifest.json with "
                     "density+rgb+weights+sample chunks) found; structural-only",
        }
        if verbose:
            print("   [functional] capture data missing - structural-only")
    else:
        try:
            report["functional"] = run_functional(
                bundle, dag_dir, out_prefix,
                psnr_threshold=psnr_threshold,
                ablation_margin=ablation_margin,
                do_ablation=do_ablation,
                gt_parity_margin=gt_parity_margin,
                check_gt_parity=check_gt_parity,
            )
            if verbose and report["functional"].get("supported"):
                f = report["functional"]
                print(f"   [functional] replay-vs-render={f['psnr_replay_vs_render']} dB (fidelity), "
                      f"replay-vs-GT={f['psnr_replay_vs_gt']} dB, "
                      f"nerfstudio-vs-GT={f['psnr_nerfstudio_vs_gt']} dB")
        except Exception as e:  # pragma: no cover - defensive
            report["functional"] = {"supported": False, "error": str(e)}
            if verbose:
                print(f"   [functional] ERROR: {e}")

    # --- Verdict ---
    reasons: List[str] = []
    structural_pass = bool(report["structural"] and report["structural"].get("pass"))
    if not structural_pass:
        reasons.append("structural checks failed")

    func = report["functional"] or {}
    functional_supported = bool(func.get("supported"))
    if functional_supported:
        rc_psnr = func.get("psnr_replay_vs_render")
        psnr_ok = rc_psnr is not None and rc_psnr >= psnr_threshold
        if not psnr_ok:
            reasons.append(f"replay-vs-render PSNR {rc_psnr} < {psnr_threshold} dB (fidelity)")
        ablations = func.get("ablations", [])
        ablation_ok = (not do_ablation) or (len(ablations) > 0 and all(a["pass"] for a in ablations))
        if do_ablation and not ablation_ok:
            failed = [a["stage"] for a in ablations if not a["pass"]]
            reasons.append(f"ablation drop < {ablation_margin} dB for: {failed or 'n/a'}")
        # GT parity: our replay must not lose quality vs nerfstudio's own render.
        gt_parity_ok = True
        if check_gt_parity and func.get("gt_parity_checked"):
            gt_parity_ok = bool(func.get("gt_parity_pass"))
            if not gt_parity_ok:
                reasons.append(
                    f"GT-parity delta {func.get('gt_parity_delta')} dB > {gt_parity_margin} dB "
                    f"(replay-vs-GT {func.get('psnr_replay_vs_gt')} vs nerfstudio-vs-GT "
                    f"{func.get('psnr_nerfstudio_vs_gt')})"
                )
        verdict_pass = structural_pass and psnr_ok and ablation_ok and gt_parity_ok
        report["verdict"] = "PASS" if verdict_pass else "FAIL"
    else:
        # Structural-only verdict when no functional capture is available.
        report["verdict"] = "PASS (structural-only)" if structural_pass else "FAIL"
        if structural_pass:
            reasons.append("functional replay skipped (no capture); structural-only verdict")

    report["verdict_reasons"] = reasons

    # --- Write report artifacts ---
    json_path = f"{out_prefix}.json"
    with open(json_path, "w") as f:
        json.dump(report, f, indent=2)
    report["report_json"] = json_path

    md_path = f"{out_prefix}.md"
    with open(md_path, "w") as f:
        f.write(_render_markdown(report))
    report["report_md"] = md_path

    if verbose:
        print(f"   [verdict] {report['verdict']}")
        print(f"   report: {json_path}")

    return report


def _render_markdown(report: Dict[str, Any]) -> str:
    lines: List[str] = ["# Operator-Graph Validation Report", ""]
    lines.append(f"**Verdict: {report['verdict']}**")
    lines.append("")
    if report.get("verdict_reasons"):
        lines.append("Notes:")
        for r in report["verdict_reasons"]:
            lines.append(f"- {r}")
        lines.append("")

    thr = report.get("thresholds", {})
    lines.append(f"- PSNR threshold (replay-vs-render, fidelity): {thr.get('psnr_threshold_db')} dB")
    lines.append(f"- Ablation margin: {thr.get('ablation_margin_db')} dB")
    if thr.get("gt_parity_margin_db") is not None:
        lines.append(f"- GT-parity margin: {thr.get('gt_parity_margin_db')} dB")
    lines.append("")

    st = report.get("structural") or {}
    lines.append("## Structural")
    lines.append(f"- Result: {'PASS' if st.get('pass') else 'FAIL'}")
    if st.get("error"):
        lines.append(f"- Error: {st['error']}")
    lines.append(f"- Nodes: {st.get('total_nodes')}  Edges: {st.get('total_edges')}")
    if st.get("stage_histogram"):
        lines.append(f"- Stage histogram: {st['stage_histogram']}")
    for c in st.get("checks", []):
        lines.append(f"  - [{'x' if c['pass'] else ' '}] {c['name']}: {c['detail']}")
    lines.append("")

    fn = report.get("functional") or {}
    lines.append("## Functional Replay")
    if not fn.get("supported"):
        lines.append(f"- Unsupported / skipped: {fn.get('error')}")
    else:
        lines.append(f"- Reference montage: `{fn.get('reference_image')}` "
                     f"(GT panel {fn.get('gt_panel_index')}, render panel "
                     f"{fn.get('render_panel_index')}, bg={fn.get('background')})")
        lines.append(f"- Rays replayed: {fn.get('n_rays')}")
        lines.append("")
        lines.append("| PSNR | dB | meaning |")
        lines.append("|---|---|---|")
        lines.append(f"| replay vs nerfstudio render | {fn.get('psnr_replay_vs_render')} | extraction fidelity (gates verdict) |")
        lines.append(f"| replay vs ground truth | {fn.get('psnr_replay_vs_gt')} | our render quality |")
        lines.append(f"| nerfstudio render vs ground truth | {fn.get('psnr_nerfstudio_vs_gt')} | nerfstudio's own quality |")
        lines.append(f"| replay vs render (checkpoint) | {fn.get('psnr_checkpoint')} | trusts captured weights |")
        lines.append("")
        if fn.get("gt_parity_checked"):
            lines.append(f"- GT parity: delta {fn.get('gt_parity_delta')} dB -> "
                         f"{'PASS' if fn.get('gt_parity_pass') else 'FAIL'}")
        lines.append(f"- Images: GT `{fn.get('ground_truth_image')}`, "
                     f"render `{fn.get('nerfstudio_render_image')}`, "
                     f"replay `{fn.get('replay_image')}`")
        lines.append("")
        lines.append("### Stage Ablations (PSNR must drop >= margin)")
        if fn.get("ablations"):
            lines.append("| stage | PSNR (dB) | drop (dB) | pass |")
            lines.append("|---|---|---|---|")
            for a in fn["ablations"]:
                lines.append(f"| {a['stage']} | {a['psnr']} | {a['psnr_drop']} | "
                             f"{'PASS' if a['pass'] else 'FAIL'} |")
        else:
            lines.append("- (ablation disabled)")
    lines.append("")
    return "\n".join(lines)


__all__ = [
    "validate_operator_graph",
    "validate_structure",
    "load_capture",
    "capture_supported",
    "run_functional",
    "PSNR_THRESHOLD_DEFAULT",
    "ABLATION_MARGIN_DEFAULT",
    "GT_PARITY_MARGIN_DEFAULT",
]