#!/usr/bin/env python3
"""Frozen Phase-3 OpProf analysis; emits aggregate, prompt-free JSON only."""

from __future__ import annotations

import argparse
import bisect
import gzip
import json
import math
import re
from collections import Counter, defaultdict
from pathlib import Path
from typing import Any

import numpy as np

import opprof_phase3_controller as common


SEED = 20260714
BOOTSTRAPS = 100_000
FAMILIES = (
    "attention",
    "moe_gemm",
    "moe_router",
    "collective",
    "sampler",
    "dense_gemm",
    "norm_elementwise",
    "kv_memory",
)
IRREGULAR_CONTROLS = (
    ("P05", "P01"),
    ("P05", "P03"),
    ("P06", "P02"),
    ("P06", "P04"),
    ("P09", "P01"),
    ("P09", "P03"),
    ("P10", "P03"),
    ("P10", "P04"),
)
AP37_MISSING_CELLS = (
    "P03/C11",
    "P05/C00",
    "P10/C00-TP2",
    "P11/C00",
)
AP37_MISSING_CONTRASTS = (
    ("P05", "P01"),
    ("P05", "P03"),
)
CAPTURE_BUCKETS = (
    1,
    2,
    4,
    8,
    16,
    24,
    32,
    40,
    48,
    56,
    64,
    72,
    80,
    88,
    96,
    104,
    112,
    120,
    128,
    136,
    144,
    152,
    160,
    168,
    176,
    184,
    192,
    200,
    208,
    216,
    224,
    232,
    240,
    248,
    256,
    272,
    288,
    304,
    320,
    336,
    352,
    368,
    384,
    400,
    416,
    432,
    448,
    464,
    480,
    496,
    512,
)


def load_json(path: Path) -> Any:
    return json.loads(path.read_text())


def numeric_sanity(values: list[float | int | None]) -> dict[str, Any]:
    finite = [float(value) for value in values if value is not None and math.isfinite(value)]
    return {
        "n": len(values),
        "finite_n": len(finite),
        "missing_n": len(values) - len(finite),
        "min": min(finite) if finite else None,
        "max": max(finite) if finite else None,
        "distinct_n": len(set(finite)),
    }


def percentile_summary(values: list[float]) -> dict[str, Any]:
    if not values:
        return {"n": 0, "mean": None, "p50": None, "p95": None, "p99": None}
    array = np.asarray(values, dtype=float)
    return {
        "n": len(values),
        "mean": float(array.mean()),
        "p50": float(np.quantile(array, 0.50)),
        "p95": float(np.quantile(array, 0.95)),
        "p99": float(np.quantile(array, 0.99)),
    }


def jsonl(path: Path) -> list[dict[str, Any]]:
    with path.open() as source:
        return [json.loads(line) for line in source]


def expected_cells() -> set[str]:
    cells = {f"P{index:02d}/C00" for index in range(1, 12)}
    for pattern in ("P01", "P03", "P06", "P10"):
        cells.update(f"{pattern}/{config}" for config in ("C10", "C01", "C11"))
    cells.add("P10/C00-TP2")
    return cells


def accepted_marker_paths(
    root: Path,
) -> tuple[list[Path], list[Path], list[Path], list[str]]:
    complete_stages = sorted(
        path.parent for path in root.glob("stages/*/stage-complete.json")
    )
    accepted_ids = []
    for stage in complete_stages:
        stage_name = stage.name
        if not (stage_name.startswith("primary-") or stage_name == "confirmations"):
            continue
        marker = load_json(stage / "stage-complete.json")
        accepted_ids.extend(str(run_id) for run_id in marker["runs"])
    if len(accepted_ids) != len(set(accepted_ids)):
        raise RuntimeError("accepted run ID appears in more than one complete stage")

    canonical = {}
    for path in root.glob("primary/*/*/run-complete.json"):
        if path.parent.name not in {"saturation", "moderate"}:
            continue
        run_id = str(load_json(path)["run_id"])
        if run_id in canonical:
            raise RuntimeError(f"duplicate canonical run marker: {run_id}")
        canonical[run_id] = path
    for path in root.glob("confirmations/*/run-complete.json"):
        run_id = str(load_json(path)["run_id"])
        if run_id in canonical:
            raise RuntimeError(f"duplicate canonical run marker: {run_id}")
        canonical[run_id] = path
    missing = sorted(set(accepted_ids) - set(canonical))
    if missing:
        raise RuntimeError(f"complete-stage run markers missing: {missing}")
    selected = [canonical[run_id] for run_id in accepted_ids]
    primary = sorted(path for path in selected if "primary" in path.parts)
    confirmations = sorted(path for path in selected if "confirmations" in path.parts)
    unaccepted = sorted(set(canonical) - set(accepted_ids))
    return primary, confirmations, complete_stages, unaccepted


def partial_verdict(has_hit: bool) -> str:
    return "PASS" if has_hit else "INCONCLUSIVE"


def run_records(run_dir: Path) -> tuple[list[dict[str, Any]], list[dict[str, Any]]]:
    result = load_json(run_dir / "client/result.json")
    t0 = int(result["t0_mono_ns"])
    start = t0 + int(60e9)
    end = t0 + int(300e9)
    stream = next((run_dir / "opprof").glob("*.jsonl"))
    all_records = []
    clean = []
    for record in jsonl(stream):
        if "step_index" not in record:
            continue
        all_records.append(record)
        if start <= int(record["submit_mono_ns"]) < end:
            clean.append(record)
    return all_records, clean


def ap36_warmup_stability(run_dir: Path) -> dict[str, Any]:
    result = load_json(run_dir / "client/result.json")
    t0 = int(result["t0_mono_ns"])
    requests = jsonl(run_dir / "client/requests.jsonl")
    completions = sum(
        bool(item["success"])
        and 0 <= float(item["completed_s"]) < float(result["warmup_seconds"])
        for item in requests
    )
    all_records, _ = run_records(run_dir)
    counts = [0, 0, 0]
    tokens = [0, 0, 0]
    for item in all_records:
        if not item["model_executed"]:
            continue
        relative_s = (int(item["submit_mono_ns"]) - t0) / 1e9
        if not 45 <= relative_s < 60:
            continue
        index = min(2, int((relative_s - 45) // 5))
        counts[index] += 1
        tokens[index] += int(item["prefill_tokens"]) + int(item["decode_tokens"])
    rates = [value / 5.0 for value in tokens]
    mean = sum(rates) / 3
    slope = (rates[2] - rates[0]) / 10.0
    drift = abs(slope) * 15 / mean if mean > 0 else math.inf
    return {
        "warmup_completions": completions,
        "step_counts": counts,
        "scheduled_tokens": tokens,
        "scheduled_token_throughput": rates,
        "mean_scheduled_token_throughput": mean,
        "slope_tokens_per_second_squared": slope,
        "normalized_drift": drift,
        "normalized_drift_limit": 0.10,
        "passes": completions >= 16
        and all(value >= 16 for value in counts)
        and drift <= 0.10,
    }


def mode(record: dict[str, Any]) -> str:
    return str(record["cudagraph"]["runtime_mode"])


def clean_block(record: dict[str, Any], t0_ns: int) -> int:
    return min(47, max(0, int((int(record["submit_mono_ns"]) - t0_ns - 60e9) // 5e9)))


def summarize_run(run_dir: Path, marker: dict[str, Any]) -> tuple[dict[str, Any], list[dict[str, Any]]]:
    client = load_json(run_dir / "client/result.json")
    requests = jsonl(run_dir / "client/requests.jsonl")
    all_records, records = run_records(run_dir)
    completed = [
        item
        for item in requests
        if item["success"] and 60 <= float(item["completed_s"]) < 300
    ]
    e2e = [float(item["completed_s"] - item["admitted_s"]) for item in completed]
    ttft = [float(item["first_token_s"] - item["admitted_s"]) for item in completed]
    tpot = [
        float(item["completed_s"] - item["first_token_s"])
        / max(1, int(item["actual_output_tokens"]) - 1)
        for item in completed
    ]
    duration_ms = [
        (int(item["complete_mono_ns"]) - int(item["submit_mono_ns"])) / 1e6
        for item in records
    ]
    useful = [int(item["prefill_tokens"]) + int(item["decode_tokens"]) for item in records]
    hit_records = [
        item
        for item in records
        if item["model_executed"]
        and item["cudagraph"]["hit"]
        and int(item["cudagraph"]["bucket_tokens"]) > 0
    ]
    model_records = [item for item in records if item["model_executed"]]
    misses = [item for item in model_records if not item["cudagraph"]["hit"]]
    eligible_misses = [
        item for item in misses if int(item["cudagraph"]["unpadded_tokens"]) <= CAPTURE_BUCKETS[-1]
    ]
    overflow = [
        item for item in misses if int(item["cudagraph"]["unpadded_tokens"]) > CAPTURE_BUCKETS[-1]
    ]
    theoretical_buckets = [
        CAPTURE_BUCKETS[bisect.bisect_left(CAPTURE_BUCKETS, int(item["cudagraph"]["unpadded_tokens"]))]
        for item in model_records
        if 0 < int(item["cudagraph"]["unpadded_tokens"]) <= CAPTURE_BUCKETS[-1]
    ]
    theoretical_unpadded = [
        int(item["cudagraph"]["unpadded_tokens"])
        for item in model_records
        if 0 < int(item["cudagraph"]["unpadded_tokens"]) <= CAPTURE_BUCKETS[-1]
    ]
    prefixes = []
    for item in records:
        for source in (item["prefix"].get("local"), item["prefix"].get("external")):
            if source:
                prefixes.append(source)
    blocks = [defaultdict(float) for _ in range(48)]
    t0 = int(client["t0_mono_ns"])
    for item, step_ms in zip(records, duration_ms, strict=True):
        block = blocks[clean_block(item, t0)]
        token = int(item["prefill_tokens"]) + int(item["decode_tokens"])
        graph = item["cudagraph"]
        block["tokens"] += token
        block["duration_ms"] += step_ms
        block["model"] += bool(item["model_executed"])
        block["miss"] += bool(item["model_executed"] and not graph["hit"])
        block["eligible_miss"] += bool(
            item["model_executed"]
            and not graph["hit"]
            and int(graph["unpadded_tokens"]) <= CAPTURE_BUCKETS[-1]
        )
        block["overflow"] += bool(
            item["model_executed"]
            and not graph["hit"]
            and int(graph["unpadded_tokens"]) > CAPTURE_BUCKETS[-1]
        )
        if item["model_executed"] and graph["hit"] and int(graph["bucket_tokens"]) > 0:
            block["padding"] += int(graph["padding_tokens"])
            block["bucket"] += int(graph["bucket_tokens"])
    denominator = sum(duration_ms)
    bucket_sum = sum(int(item["cudagraph"]["bucket_tokens"]) for item in hit_records)
    pad_sum = sum(int(item["cudagraph"]["padding_tokens"]) for item in hit_records)
    model_n = len(model_records)
    modes = Counter(mode(item) for item in records)
    completion_times = [
        float(item["completed_s"])
        for item in requests
        if item["success"]
    ]
    clean_completion_rates = [
        sum(start <= value < start + 10 for value in completion_times) / 10
        for start in np.arange(220, 300, 10)
    ]
    clean_waiting = [
        item["queues"]["waiting"]
        for item in records
        if t0 + int(220e9) <= int(item["submit_mono_ns"]) < t0 + int(300e9)
    ]
    clean_completion_median = float(np.median(clean_completion_rates))
    clean_waiting_median = float(np.median(clean_waiting)) if clean_waiting else None
    recoveries = []
    profiles = client.get("profiles", [])
    for index, profile in enumerate(profiles):
        recovery_end = (
            float(profiles[index + 1]["start_call_s"])
            if index + 1 < len(profiles)
            else float(client["admission_stop_s"])
        )
        tail_start = recovery_end - 10
        completion_rate = (
            sum(tail_start <= value < recovery_end for value in completion_times) / 10
        )
        tail_waiting = [
            item["queues"]["waiting"]
            for item in all_records
            if t0 + int(tail_start * 1e9)
            <= int(item["submit_mono_ns"])
            < t0 + int(recovery_end * 1e9)
        ]
        waiting_median = float(np.median(tail_waiting)) if tail_waiting else None
        rate_ok = (
            completion_rate == clean_completion_median
            if clean_completion_median == 0
            else abs(completion_rate / clean_completion_median - 1) <= 0.10
        )
        waiting_ok = (
            waiting_median == clean_waiting_median
            if clean_waiting_median == 0
            else (
                waiting_median is not None
                and clean_waiting_median is not None
                and abs(waiting_median / clean_waiting_median - 1) <= 0.10
            )
        )
        recoveries.append(
            {
                "window": index + 1,
                "tail_start_s": tail_start,
                "tail_end_s": recovery_end,
                "completion_rate_rps": completion_rate,
                "clean_c_completion_rate_median_rps": clean_completion_median,
                "waiting_median": waiting_median,
                "clean_c_waiting_median": clean_waiting_median,
                "completion_rate_within_10pct": rate_ok,
                "waiting_within_10pct": waiting_ok,
                "valid": rate_ok and waiting_ok,
            }
        )
    summary = {
        "run_id": marker["run_id"],
        "pattern": marker["pattern"],
        "config": marker["config"],
        "load": client["load_point"],
        "drain_seconds": float(client["drain_seconds"]),
        "drain_quarantined": bool(marker["drain_quarantined"]),
        "clean": client["clean"],
        "requests_completed": len(completed),
        "latency_s": {
            "e2e": percentile_summary(e2e),
            "ttft": percentile_summary(ttft),
            "tpot": percentile_summary(tpot),
        },
        "layer1": {
            "records_all": len(all_records),
            "records_clean": len(records),
            "step_duration_ms": percentile_summary(duration_ms),
            "scheduled_tokens_per_step": percentile_summary([float(value) for value in useful]),
            "requests_per_step": percentile_summary(
                [float(item["scheduled_requests"]) for item in records]
            ),
            "prefill_tokens": sum(int(item["prefill_tokens"]) for item in records),
            "decode_tokens": sum(int(item["decode_tokens"]) for item in records),
            "token_efficiency_per_ms": sum(useful) / denominator if denominator else None,
            "queue_waiting_mean": float(np.mean([item["queues"]["waiting"] for item in records])),
            "queue_waiting_max": max((item["queues"]["waiting"] for item in records), default=0),
            "kv_usage_mean": float(np.mean([item["kv"]["usage"] for item in records])),
            "kv_usage_max": max((item["kv"]["usage"] for item in records), default=0),
            "preemptions": sum(int(item["preemptions"]) for item in records),
            "prefix_query_hit_ratio": (
                sum(int(item["hits"]) for item in prefixes)
                / sum(int(item["queries"]) for item in prefixes)
                if sum(int(item["queries"]) for item in prefixes)
                else 0.0
            ),
            "mode_counts": dict(sorted(modes.items())),
            "mode_shares": {key: value / len(records) for key, value in sorted(modes.items())},
        },
        "waste": {
            "padding_fraction": pad_sum / bucket_sum if bucket_sum else 0.0,
            "padded_tokens_per_useful_token": pad_sum / sum(useful) if sum(useful) else 0.0,
            "graph_miss_rate": len(misses) / model_n if model_n else 0.0,
            "eligible_miss_rate": len(eligible_misses) / model_n if model_n else 0.0,
            "overflow_rate": len(overflow) / model_n if model_n else 0.0,
            "bucket_slack": (
                (sum(theoretical_buckets) - sum(theoretical_unpadded))
                / sum(theoretical_buckets)
                if theoretical_buckets
                else 0.0
            ),
            "mixed_interference": None,
            "moe_layer_cv": None,
        },
        "trace_files": len(marker.get("traces", [])),
        "profile_recovery": recoveries,
        "profile_recovery_valid": all(item["valid"] for item in recoveries),
        "layer2_missing_after_controller_cleanup": bool(
            marker.get("layer2_missing_after_controller_cleanup")
        ),
        "blocks": [dict(block) for block in blocks],
    }
    return summary, records


def manifest_raggedness(path: Path, cohort: int) -> tuple[float, list[tuple[float, float]]]:
    lengths = []
    with path.open() as source:
        for line in source:
            lengths.append(int(json.loads(line)["input_tokens"]))
    pieces = []
    for start in range(0, len(lengths) - cohort + 1, cohort):
        group = lengths[start : start + cohort]
        denominator = float(cohort * max(group))
        pieces.append((denominator - sum(group), denominator))
    numerator = sum(item[0] for item in pieces)
    denominator = sum(item[1] for item in pieces)
    return numerator / denominator, pieces


EXECUTE = re.compile(r"execute_context_\d+\((\d+)\)_generation_\d+\((\d+)\)")


def trace_steps(path: Path, layer1: list[dict[str, Any]]) -> dict[str, Any]:
    opener = gzip.open if path.suffix == ".gz" else open
    with opener(path, "rt", encoding="utf-8") as source:
        payload = json.load(source)
    base = int(payload["baseTimeNanoseconds"])
    events = payload["traceEvents"]
    cpu = []
    gpu_by_external: dict[int, list[dict[str, Any]]] = defaultdict(list)
    kernels = []
    for event in events:
        name = str(event.get("name", ""))
        if event.get("cat") == "user_annotation" and EXECUTE.fullmatch(name):
            cpu.append(event)
        elif event.get("cat") == "gpu_user_annotation" and EXECUTE.fullmatch(name):
            gpu_by_external[int(event.get("args", {}).get("External id", -1))].append(event)
        elif event.get("cat") == "kernel":
            kernels.append(event)
    cpu.sort(key=lambda item: float(item["ts"]))
    if len(cpu) != 8:
        raise RuntimeError(f"active execute count {len(cpu)} != 8: {path}")
    walls = [int(item["submit_wall_ns"]) for item in layer1]
    joined = []
    used_steps = set()
    unmatched: dict[str, float] = defaultdict(float)
    for event in cpu:
        external = int(event["args"]["External id"])
        candidates = gpu_by_external.get(external, [])
        if not candidates:
            raise RuntimeError(f"GPU execute annotation absent: {path}: {external}")
        gpu = max(candidates, key=lambda item: float(item.get("dur", 0)))
        start = float(gpu["ts"])
        end = start + float(gpu["dur"])
        durations: dict[str, float] = defaultdict(float)
        for kernel in kernels:
            midpoint = float(kernel["ts"]) + float(kernel.get("dur", 0)) / 2
            if start <= midpoint <= end:
                kernel_name = str(kernel.get("name", ""))
                duration = float(kernel.get("dur", 0))
                family = common.classify_kernel(kernel_name)
                durations[family] += duration
                if family == "other":
                    unmatched[kernel_name] += duration
        total = sum(durations.values())
        if total <= 0:
            raise RuntimeError(f"active execute has no kernels: {path}: {external}")
        wall = base + int(float(event["ts"]) * 1000)
        index = bisect.bisect_left(walls, wall)
        choices = layer1[max(0, index - 3) : index + 3]
        record = min(choices, key=lambda item: abs(int(item["submit_wall_ns"]) - wall))
        match = EXECUTE.fullmatch(str(event["name"]))
        assert match
        expected = (int(match.group(1)), int(match.group(2)))
        actual = (int(record["prefill_tokens"]), int(record["decode_tokens"]))
        delta_ms = (int(record["submit_wall_ns"]) - wall) / 1e6
        if actual != expected or abs(delta_ms) > 100 or record["step_index"] in used_steps:
            raise RuntimeError(
                f"ambiguous Layer-1 join: {path}: expected={expected} actual={actual} "
                f"delta_ms={delta_ms}"
            )
        used_steps.add(record["step_index"])
        shares = {family: durations.get(family, 0.0) / total for family in FAMILIES}
        shares["other"] = durations.get("other", 0.0) / total
        joined.append(
            {
                "step_index": int(record["step_index"]),
                "join_delta_ms": delta_ms,
                "prefill_tokens": actual[0],
                "decode_tokens": actual[1],
                "scheduled_requests": int(record["scheduled_requests"]),
                "decode_batch_size": int(record["decode_batch_size"]),
                "runtime_mode": mode(record),
                "duration_us": dict(durations),
                "shares": shares,
                "classifiable_fraction": 1.0 - shares["other"],
            }
        )
    aggregate = defaultdict(float)
    for step in joined:
        for family, duration in step["duration_us"].items():
            aggregate[family] += duration
    total = sum(aggregate.values())
    attention_subduration = defaultdict(float)
    mode_duration: dict[str, dict[str, float]] = defaultdict(lambda: defaultdict(float))
    mode_steps = Counter()
    for step in joined:
        if step["prefill_tokens"] and step["decode_tokens"]:
            attention_key = "attention_mixed"
        elif step["prefill_tokens"]:
            attention_key = "attention_prefill"
        else:
            attention_key = "attention_decode"
        attention_subduration[attention_key] += step["duration_us"].get("attention", 0.0)
        mode_steps[step["runtime_mode"]] += 1
        for family, duration in step["duration_us"].items():
            mode_duration[step["runtime_mode"]][family] += duration
    return {
        "path": path.name,
        "steps": joined,
        "shares": {family: aggregate.get(family, 0.0) / total for family in FAMILIES},
        "other_share": aggregate.get("other", 0.0) / total,
        "classifiable_fraction": 1.0 - aggregate.get("other", 0.0) / total,
        "attention_subshares": {
            key: attention_subduration.get(key, 0.0) / total
            for key in ("attention_prefill", "attention_decode", "attention_mixed")
        },
        "mode_steps": dict(sorted(mode_steps.items())),
        "mode_shares": {
            key: {
                family: durations.get(family, 0.0) / sum(durations.values())
                for family in FAMILIES
            }
            for key, durations in sorted(mode_duration.items())
        },
        "top_unmatched": [
            {"name": name, "duration_us": duration}
            for name, duration in sorted(unmatched.items(), key=lambda item: item[1], reverse=True)[
                :20
            ]
        ],
    }


def smd(profile: list[float], clean: list[float]) -> float:
    a = np.asarray(profile, dtype=float)
    b = np.asarray(clean, dtype=float)
    if np.all(a == a[0]) and np.all(b == b[0]):
        return 0.0 if a[0] == b[0] else math.inf
    denominator = math.sqrt((float(a.var(ddof=1)) + float(b.var(ddof=1))) / 2)
    if denominator == 0:
        return math.inf
    return float((a.mean() - b.mean()) / denominator)


def representativeness(window: dict[str, Any], clean_c: list[dict[str, Any]]) -> dict[str, Any]:
    steps = window["steps"]
    features: dict[str, tuple[list[float], list[float]]] = {
        "scheduled_tokens": (
            [step["prefill_tokens"] + step["decode_tokens"] for step in steps],
            [item["prefill_tokens"] + item["decode_tokens"] for item in clean_c],
        ),
        "prefill_fraction": (
            [
                step["prefill_tokens"]
                / max(1, step["prefill_tokens"] + step["decode_tokens"])
                for step in steps
            ],
            [
                item["prefill_tokens"]
                / max(1, item["prefill_tokens"] + item["decode_tokens"])
                for item in clean_c
            ],
        ),
        "decode_batch_size": (
            [step["decode_batch_size"] for step in steps],
            [item["decode_batch_size"] for item in clean_c],
        ),
    }
    modes = ("FULL", "PIECEWISE", "NONE")
    for key in modes:
        features[f"mode_{key}"] = (
            [float(step["runtime_mode"] == key) for step in steps],
            [float(mode(item) == key) for item in clean_c],
        )
    values = {key: smd(*pair) for key, pair in features.items()}
    return {
        "smd": values,
        "valid": all(math.isfinite(value) and abs(value) <= 0.25 for value in values.values()),
    }


def convex_hull(points: list[tuple[float, float]]) -> list[tuple[float, float]]:
    unique = sorted(set(points))
    if len(unique) <= 1:
        return unique

    def cross(origin, left, right):
        return (left[0] - origin[0]) * (right[1] - origin[1]) - (
            left[1] - origin[1]
        ) * (right[0] - origin[0])

    lower = []
    for point in unique:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], point) <= 0:
            lower.pop()
        lower.append(point)
    upper = []
    for point in reversed(unique):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], point) <= 0:
            upper.pop()
        upper.append(point)
    return lower[:-1] + upper[:-1]


def inside_convex(hull: list[tuple[float, float]], point: tuple[float, float]) -> bool:
    if not hull:
        return False
    if len(hull) == 1:
        return point == hull[0]
    if len(hull) == 2:
        left, right = hull
        cross = (right[0] - left[0]) * (point[1] - left[1]) - (
            right[1] - left[1]
        ) * (point[0] - left[0])
        return (
            abs(cross) <= 1e-9
            and min(left[0], right[0]) <= point[0] <= max(left[0], right[0])
            and min(left[1], right[1]) <= point[1] <= max(left[1], right[1])
        )
    signs = []
    for index, left in enumerate(hull):
        right = hull[(index + 1) % len(hull)]
        cross = (right[0] - left[0]) * (point[1] - left[1]) - (
            right[1] - left[1]
        ) * (point[0] - left[0])
        if abs(cross) > 1e-9:
            signs.append(math.copysign(1, cross))
    return not signs or min(signs) == max(signs)


def fit_nonnegative_robust(rows: list[tuple[float, float, float]]):
    values = np.asarray(rows, dtype=float)
    x = np.log1p(values[:, 0])
    n = np.log1p(values[:, 1])
    y = values[:, 2]
    x_knots = np.quantile(x, [0.25, 0.5, 0.75])
    n_knots = np.quantile(n, [0.25, 0.5, 0.75])

    def design(raw_x, raw_n):
        a = np.log1p(np.asarray(raw_x, dtype=float))
        b = np.log1p(np.asarray(raw_n, dtype=float))
        columns = [np.ones_like(a), a, b, a * b]
        columns.extend(np.maximum(0, a - knot) for knot in x_knots)
        columns.extend(np.maximum(0, b - knot) for knot in n_knots)
        return np.column_stack(columns)

    matrix = design(values[:, 0], values[:, 1])
    weights = np.ones(len(y))
    coef = np.zeros(matrix.shape[1])
    for _ in range(20):
        weighted = np.sqrt(weights)
        coef = np.linalg.lstsq(matrix * weighted[:, None], y * weighted, rcond=None)[0]
        coef = np.maximum(0, coef)
        residual = y - matrix @ coef
        scale = 1.4826 * np.median(np.abs(residual - np.median(residual))) + 1e-9
        weights = np.minimum(1.0, 1.345 * scale / np.maximum(np.abs(residual), 1e-9))

    def predict(raw_x: float, raw_n: float) -> float:
        return float(max(0.0, design([raw_x], [raw_n])[0] @ coef))

    return predict, convex_hull(
        [(float(row[0]), float(row[1])) for row in values]
    )


def add_mixed_interference(
    summaries: dict[str, dict[str, Any]],
    records: dict[str, list[dict[str, Any]]],
) -> None:
    groups: dict[tuple[str, str], list[str]] = defaultdict(list)
    for run_id, summary in summaries.items():
        if run_id.endswith("-confirmation"):
            continue
        groups[(summary["config"], summary["load"])].append(run_id)
    for run_ids in groups.values():
        for target in run_ids:
            pattern = summaries[target]["pattern"]
            training = [item for item in run_ids if summaries[item]["pattern"] != pattern]
            prefill = []
            decode = []
            zero = []
            for run_id in training:
                for item in records[run_id]:
                    p = int(item["prefill_tokens"])
                    d = int(item["decode_tokens"])
                    n = int(item["scheduled_requests"])
                    duration = (item["complete_mono_ns"] - item["submit_mono_ns"]) / 1e6
                    if p > 0 and d == 0:
                        prefill.append((p, n, duration))
                    elif d > 0 and p == 0:
                        decode.append((d, n, duration))
                    elif p == 0 and d == 0:
                        zero.append(duration)
            if len(prefill) < 30 or len(decode) < 30:
                continue
            fp, p_support = fit_nonnegative_robust(prefill)
            fd, d_support = fit_nonnegative_robust(decode)
            alpha = float(np.median(zero)) if zero else 0.0
            supported = []
            for item in records[target]:
                p = int(item["prefill_tokens"])
                d = int(item["decode_tokens"])
                n = int(item["scheduled_requests"])
                if not (p > 0 and d > 0):
                    continue
                if not inside_convex(p_support, (p, n)) or not inside_convex(
                    d_support, (d, n)
                ):
                    continue
                predicted = fp(p, n) + fd(d, n) - alpha
                if predicted <= 0:
                    continue
                observed = (item["complete_mono_ns"] - item["submit_mono_ns"]) / 1e6
                supported.append((item, observed - predicted, predicted))
            summaries[target]["waste"]["mixed_supported_steps"] = len(supported)
            summaries[target]["waste"]["mixed_total_steps"] = sum(
                item["prefill_tokens"] > 0 and item["decode_tokens"] > 0
                for item in records[target]
            )
            if len(supported) < 30:
                continue
            numerator = sum(item[1] for item in supported)
            denominator = sum(item[2] for item in supported)
            summaries[target]["waste"]["mixed_interference"] = numerator / denominator
            blocks = summaries[target]["blocks"]
            t0 = int(load_json(Path(summaries[target]["run_dir"]) / "client/result.json")["t0_mono_ns"])
            for record, residual, predicted in supported:
                block = blocks[clean_block(record, t0)]
                block["mix_residual"] = block.get("mix_residual", 0.0) + residual
                block["mix_predicted"] = block.get("mix_predicted", 0.0) + predicted


def ratio_blocks(summary: dict[str, Any], metric: str) -> np.ndarray | None:
    keys = {
        "padding_fraction": ("padding", "bucket"),
        "graph_miss_rate": ("miss", "model"),
        "overflow_rate": ("overflow", "model"),
        "mixed_interference": ("mix_residual", "mix_predicted"),
        "efficiency": ("tokens", "duration_ms"),
    }
    numerator, denominator = keys[metric]
    values = np.asarray(
        [[block.get(numerator, 0.0), block.get(denominator, 0.0)] for block in summary["blocks"]],
        dtype=float,
    )
    if values[:, 1].sum() <= 0:
        return None
    return values


def bootstrap_difference(
    left: np.ndarray, right: np.ndarray, rng: np.random.Generator
) -> dict[str, float]:
    point = left[:, 0].sum() / left[:, 1].sum() - right[:, 0].sum() / right[:, 1].sum()
    draws = np.empty(BOOTSTRAPS)
    chunk = 5000
    for start in range(0, BOOTSTRAPS, chunk):
        count = min(chunk, BOOTSTRAPS - start)
        li = rng.integers(0, len(left), size=(count, len(left)))
        ri = rng.integers(0, len(right), size=(count, len(right)))
        ln = left[li, 0].sum(axis=1)
        ld = left[li, 1].sum(axis=1)
        rn = right[ri, 0].sum(axis=1)
        rd = right[ri, 1].sum(axis=1)
        draws[start : start + count] = ln / np.maximum(ld, 1e-12) - rn / np.maximum(rd, 1e-12)
    p = min(1.0, 2 * min(float(np.mean(draws <= 0)), float(np.mean(draws >= 0))))
    return {
        "point": float(point),
        "ci95_low": float(np.quantile(draws, 0.025)),
        "ci95_high": float(np.quantile(draws, 0.975)),
        "simultaneous_low": float(np.quantile(draws, 0.05 / 16)),
        "simultaneous_high": float(np.quantile(draws, 1 - 0.05 / 16)),
        "p": p,
    }


def holm(results: list[dict[str, Any]], total_tests: int | None = None) -> None:
    ordered = sorted(results, key=lambda item: item["p"])
    count = total_tests or len(ordered)
    running = 0.0
    for index, item in enumerate(ordered):
        running = max(running, min(1.0, (count - index) * item["p"]))
        item["p_holm"] = running


def permutation_p(
    left: list[dict[str, Any]],
    right: list[dict[str, Any]],
    family_a: str,
    family_b: str,
    orientation: float,
    rng: np.random.Generator,
) -> float:
    observed_windows = []
    null_windows = []
    for window in range(2):
        a = np.asarray(
            [
                orientation * (step["shares"][family_a] - step["shares"][family_b])
                for step in left[window]["steps"]
            ]
        )
        b = np.asarray(
            [
                orientation * (step["shares"][family_a] - step["shares"][family_b])
                for step in right[window]["steps"]
            ]
        )
        observed_windows.append(float(a.mean() - b.mean()))
        pooled = np.concatenate([a, b])
        values = np.empty(BOOTSTRAPS)
        for start in range(0, BOOTSTRAPS, 5000):
            count = min(5000, BOOTSTRAPS - start)
            random = rng.random((count, 16))
            indices = np.argpartition(random, 8, axis=1)[:, :8]
            selected = np.take_along_axis(np.broadcast_to(pooled, (count, 16)), indices, axis=1)
            values[start : start + count] = 2 * selected.mean(axis=1) - pooled.mean() * 2
        null_windows.append(values)
    observed = min(observed_windows)
    null = np.minimum(null_windows[0], null_windows[1])
    return float((np.count_nonzero(null >= observed) + 1) / (BOOTSTRAPS + 1))


def kendall_tau_b(left: list[float], right: list[float]) -> float:
    concordant = discordant = tie_left = tie_right = 0
    for i in range(len(left)):
        for j in range(i + 1, len(left)):
            a = np.sign(left[i] - left[j])
            b = np.sign(right[i] - right[j])
            if a == 0 and b != 0:
                tie_left += 1
            elif b == 0 and a != 0:
                tie_right += 1
            elif a * b > 0:
                concordant += 1
            elif a * b < 0:
                discordant += 1
    denominator = math.sqrt(
        (concordant + discordant + tie_left) * (concordant + discordant + tie_right)
    )
    return (concordant - discordant) / denominator if denominator else 0.0


def ranked_families(shares: dict[str, float]) -> list[dict[str, Any]]:
    ordered = sorted(FAMILIES, key=lambda family: shares[family], reverse=True)
    result = []
    rank = 1
    group_top = None
    for index, family in enumerate(ordered):
        if group_top is None or group_top - shares[family] > 0.01:
            rank = index + 1
            group_top = shares[family]
        result.append({"family": family, "share": shares[family], "rank": rank})
    return result


def operator_window_valid(
    window: dict[str, Any], summary: dict[str, Any]
) -> bool:
    representative = window.get("representativeness")
    return bool(
        window["classifiable_fraction"] >= 0.70
        and representative
        and representative["valid"]
        and summary["profile_recovery_valid"]
    )


def analyze(root: Path, private: Path) -> dict[str, Any]:
    markers, confirmation_markers, complete_stages, unaccepted_markers = (
        accepted_marker_paths(root)
    )
    if len(markers) != 40 or confirmation_markers:
        raise RuntimeError(
            f"A-P3-7 run count before analysis: primary={len(markers)} "
            f"confirm={len(confirmation_markers)}"
        )
    summaries: dict[str, dict[str, Any]] = {}
    layer_records: dict[str, list[dict[str, Any]]] = {}
    for path in markers + confirmation_markers:
        marker = load_json(path)
        run_id = marker["run_id"]
        summary, records = summarize_run(path.parent, marker)
        summary["run_dir"] = str(path.parent)
        summaries[run_id] = summary
        layer_records[run_id] = records
    cell_loads: dict[str, set[str]] = defaultdict(set)
    for summary in summaries.values():
        cell_loads[f"{summary['pattern']}/{summary['config']}"].add(summary["load"])
    complete_cells = sorted(
        cell for cell, loads in cell_loads.items() if loads == {"saturation", "moderate"}
    )
    missing_cells = sorted(expected_cells() - set(complete_cells))
    if missing_cells != sorted(AP37_MISSING_CELLS):
        raise RuntimeError(f"A-P3-7 missing-cell mismatch: {missing_cells}")
    completed_c00_patterns = sorted(
        cell.split("/", 1)[0] for cell in complete_cells if cell.endswith("/C00")
    )
    ragged = {}
    ragged_pieces = {}
    for pattern in (f"P{index:02d}" for index in range(1, 12)):
        values = {}
        for cohort in (32, 64, 128):
            value, pieces = manifest_raggedness(private / f"{pattern}.jsonl", cohort)
            values[f"R{cohort}"] = value
            if cohort == 64:
                ragged_pieces[pattern] = np.asarray(pieces)
        ragged[pattern] = values
    for summary in summaries.values():
        summary["waste"].update(ragged[summary["pattern"]])
    add_mixed_interference(summaries, layer_records)

    layer2: dict[str, list[dict[str, Any]]] = {}
    layer2_issues = []
    for run_id, summary in summaries.items():
        run_dir = Path(summary["run_dir"])
        traces = sorted((run_dir / "traces").glob("*.pt.trace.json*"))
        if not traces:
            if not summary["layer2_missing_after_controller_cleanup"] and not run_id.endswith("confirmation"):
                layer2_issues.append(f"unexpected missing traces: {run_id}")
            continue
        all_layer, _ = run_records(run_dir)
        windows = [trace_steps(path, all_layer) for path in traces]
        if summary["config"] != "C00-TP2":
            clean_result = load_json(run_dir / "client/result.json")
            t0 = int(clean_result["t0_mono_ns"])
            clean_c = [
                item
                for item in all_layer
                if t0 + int(220e9) <= int(item["submit_mono_ns"]) < t0 + int(300e9)
            ]
            for window in windows:
                window["representativeness"] = representativeness(window, clean_c)
        layer2[run_id] = windows

    all_operator_windows = {
        run_id: [
            {
                "trace": window["path"],
                "shares": window["shares"],
                "other_share": window["other_share"],
                "classifiable_fraction": window["classifiable_fraction"],
                "attention_subshares": window["attention_subshares"],
                "mode_steps": window["mode_steps"],
                "mode_shares": window["mode_shares"],
                "top_unmatched": window["top_unmatched"],
                "representativeness": window.get("representativeness"),
            }
            for window in windows
        ]
        for run_id, windows in sorted(layer2.items())
    }

    patterns = [f"P{i:02d}" for i in range(1, 12)]
    primary_ids = [f"{pattern}-C00-moderate" for pattern in completed_c00_patterns]
    operator_table = {}
    operator_share_table = {}
    operator_mode_segments = {}
    invalid_primary_windows = []
    for pattern in patterns:
        operator_share_table[pattern] = {}
        for load in ("saturation", "moderate"):
            run_id = f"{pattern}-C00-{load}"
            windows = layer2.get(run_id, [])
            rows = []
            for index, window in enumerate(windows, 1):
                valid = operator_window_valid(window, summaries[run_id])
                rows.append(
                    {
                        "window": index,
                        "shares": window["shares"],
                        "ranking": ranked_families(window["shares"]),
                        "classifiable_fraction": window["classifiable_fraction"],
                        "representativeness": window["representativeness"],
                        "attention_subshares": window["attention_subshares"],
                        "mode_steps": window["mode_steps"],
                        "mode_shares": window["mode_shares"],
                        "top_unmatched": window["top_unmatched"],
                        "valid": valid,
                    }
                )
            mean_shares = (
                {
                    family: float(
                        np.mean([window["shares"][family] for window in windows])
                    )
                    for family in FAMILIES
                }
                if windows
                else None
            )
            status = (
                "EVALUABLE"
                if len(rows) == 2 and all(row["valid"] for row in rows)
                else "NOT_EVALUABLE"
            )
            reason = None
            if run_id not in summaries:
                reason = "pattern/config cell missing under A-P3-7"
            elif len(rows) != 2:
                reason = f"accepted trace windows {len(rows)}/2"
            elif status != "EVALUABLE":
                reason = "classifiability, representativeness, or recovery gate failed"
            operator_share_table[pattern][load] = {
                "run_id": run_id,
                "status": status,
                "reason": reason,
                "window_count": len(rows),
                "valid_window_count": sum(row["valid"] for row in rows),
                "mean_shares": mean_shares,
                "mean_ranking": ranked_families(mean_shares) if mean_shares else None,
                "windows": rows,
            }
            if run_id in primary_ids:
                if len(rows) != 2:
                    invalid_primary_windows.append(
                        f"{run_id}: trace count {len(rows)}"
                    )
                else:
                    operator_table[run_id] = rows
                    invalid_primary_windows.extend(
                        f"{run_id}: window {row['window']}"
                        for row in rows
                        if not row["valid"]
                    )
            if not windows:
                continue
            by_mode: dict[str, list[dict[str, Any]]] = defaultdict(list)
            for window in windows:
                for step in window["steps"]:
                    by_mode[step["runtime_mode"]].append(step)
            operator_mode_segments[run_id] = {}
            for runtime_mode, steps in sorted(by_mode.items()):
                if len(steps) < 8:
                    continue
                durations = defaultdict(float)
                for step in steps:
                    for family, duration in step["duration_us"].items():
                        durations[family] += duration
                total = sum(durations.values())
                operator_mode_segments[run_id][runtime_mode] = {
                    "steps": len(steps),
                    "shares": {
                        family: durations.get(family, 0.0) / total
                        for family in FAMILIES
                    },
                    "other_share": durations.get("other", 0.0) / total,
                }

    rng = np.random.default_rng(SEED)
    inversions = []
    total_ranking_tests = math.comb(11, 2) * math.comb(len(FAMILIES), 2)
    for pi, left_pattern in enumerate(patterns):
        left_id = f"{left_pattern}-C00-moderate"
        if len(operator_table.get(left_id, [])) != 2 or any(
            not item["valid"] for item in operator_table.get(left_id, [])
        ):
            continue
        for right_pattern in patterns[pi + 1 :]:
            right_id = f"{right_pattern}-C00-moderate"
            if len(operator_table.get(right_id, [])) != 2 or any(
                not item["valid"] for item in operator_table.get(right_id, [])
            ):
                continue
            for ai, family_a in enumerate(FAMILIES):
                for family_b in FAMILIES[ai + 1 :]:
                    left_gaps = [
                        window["shares"][family_a] - window["shares"][family_b]
                        for window in layer2[left_id]
                    ]
                    right_gaps = [
                        window["shares"][family_a] - window["shares"][family_b]
                        for window in layer2[right_id]
                    ]
                    orientation = None
                    if min(left_gaps) >= 0.05 and max(right_gaps) <= -0.05:
                        orientation = 1.0
                    elif max(left_gaps) <= -0.05 and min(right_gaps) >= 0.05:
                        orientation = -1.0
                    if orientation is None:
                        continue
                    p = permutation_p(
                        layer2[left_id],
                        layer2[right_id],
                        family_a,
                        family_b,
                        orientation,
                        rng,
                    )
                    inversions.append(
                        {
                            "left": left_pattern,
                            "right": right_pattern,
                            "family_a": family_a,
                            "family_b": family_b,
                            "left_gaps": left_gaps,
                            "right_gaps": right_gaps,
                            "orientation": orientation,
                            "p": p,
                        }
                    )
    holm(inversions, total_tests=total_ranking_tests)
    accepted_inversions = [item for item in inversions if item["p_holm"] < 0.05]

    tau = []
    for index, left in enumerate(patterns):
        left_id = f"{left}-C00-moderate"
        if len(operator_table.get(left_id, [])) != 2 or any(
            not item["valid"] for item in operator_table.get(left_id, [])
        ):
            continue
        left_mean = {
            family: float(np.mean([window["shares"][family] for window in layer2[left_id]]))
            for family in FAMILIES
        }
        left_rank = {item["family"]: item["rank"] for item in ranked_families(left_mean)}
        for right in patterns[index + 1 :]:
            right_id = f"{right}-C00-moderate"
            if len(operator_table.get(right_id, [])) != 2 or any(
                not item["valid"] for item in operator_table.get(right_id, [])
            ):
                continue
            right_mean = {
                family: float(
                    np.mean([window["shares"][family] for window in layer2[right_id]])
                )
                for family in FAMILIES
            }
            right_rank = {
                item["family"]: item["rank"] for item in ranked_families(right_mean)
            }
            tau.append(
                {
                    "left": left,
                    "right": right,
                    "tau_b": kendall_tau_b(
                        [left_rank[family] for family in FAMILIES],
                        [right_rank[family] for family in FAMILIES],
                    ),
                }
            )

    waste_contrasts = []
    contrast_status = []
    by_metric: dict[str, list[dict[str, Any]]] = defaultdict(list)
    for irregular, control in IRREGULAR_CONTROLS:
        left_id = f"{irregular}-C00-moderate"
        right_id = f"{control}-C00-moderate"
        missing = [
            run_id.removesuffix("-moderate")
            for run_id in (left_id, right_id)
            if run_id not in summaries
        ]
        if missing:
            contrast_status.append(
                {
                    "irregular": irregular,
                    "control": control,
                    "evaluable": False,
                    "verdict": "NOT_EVALUABLE",
                    "missing_cells": missing,
                    "efficiency_loss": None,
                    "metrics": [],
                    "passing_metrics": [],
                    "passes_h1b": False,
                }
            )
            continue
        left = summaries[left_id]
        right = summaries[right_id]
        efficiency_loss = 1 - (
            left["layer1"]["token_efficiency_per_ms"]
            / right["layer1"]["token_efficiency_per_ms"]
        )
        metric_results = []
        for metric in ("padding_fraction", "graph_miss_rate", "overflow_rate", "mixed_interference"):
            a = ratio_blocks(left, metric)
            b = ratio_blocks(right, metric)
            if a is None or b is None:
                result = {
                    "point": None,
                    "ci95_low": None,
                    "ci95_high": None,
                    "simultaneous_low": None,
                    "simultaneous_high": None,
                    "p": 1.0,
                    "available": False,
                }
            else:
                result = bootstrap_difference(a, b, rng)
                result["available"] = True
            result.update(
                {
                    "irregular": irregular,
                    "control": control,
                    "metric": metric,
                    "efficiency_loss": efficiency_loss,
                    "evaluable": True,
                }
            )
            by_metric[metric].append(result)
            waste_contrasts.append(result)
            metric_results.append(result)
        result = bootstrap_difference(ragged_pieces[irregular], ragged_pieces[control], rng)
        result["available"] = True
        result.update(
            {
                "irregular": irregular,
                "control": control,
                "metric": "R64",
                "efficiency_loss": efficiency_loss,
                "evaluable": True,
            }
        )
        by_metric["R64"].append(result)
        waste_contrasts.append(result)
        metric_results.append(result)
        result = {
            "irregular": irregular,
            "control": control,
            "metric": "moe_layer_cv",
            "efficiency_loss": efficiency_loss,
            "evaluable": True,
            "point": None,
            "ci95_low": None,
            "ci95_high": None,
            "simultaneous_low": None,
            "simultaneous_high": None,
            "p": 1.0,
            "available": False,
        }
        by_metric["moe_layer_cv"].append(result)
        waste_contrasts.append(result)
        metric_results.append(result)
        contrast_status.append(
            {
                "irregular": irregular,
                "control": control,
                "evaluable": True,
                "verdict": None,
                "missing_cells": [],
                "efficiency_loss": efficiency_loss,
                "metrics": metric_results,
                "passing_metrics": [],
                "passes_h1b": False,
            }
        )
    for values in by_metric.values():
        holm(values, total_tests=len(IRREGULAR_CONTROLS))
    thresholds = {
        "padding_fraction": 0.05,
        "graph_miss_rate": 0.10,
        "overflow_rate": 0.10,
        "R64": 0.15,
        "mixed_interference": 0.10,
        "moe_layer_cv": 0.15,
    }
    for item in waste_contrasts:
        item["material"] = (
            item["point"] is not None
            and item["point"] >= thresholds[item["metric"]]
            and item["simultaneous_low"] is not None
            and item["simultaneous_low"] > 0
            and item["p_holm"] < 0.05
        )
        item["coincident_efficiency_or_residual"] = (
            item["efficiency_loss"] >= 0.05
            or (
                summaries[f"{item['irregular']}-C00-moderate"]["waste"].get(
                    "mixed_interference"
                )
                or -math.inf
            )
            > 0
        )
        item["passes_h1b"] = item["material"] and item["coincident_efficiency_or_residual"]
    for contrast in contrast_status:
        if not contrast["evaluable"]:
            continue
        hits = [item for item in contrast["metrics"] if item["passes_h1b"]]
        contrast["passing_metrics"] = [item["metric"] for item in hits]
        contrast["passes_h1b"] = bool(hits)
        contrast["verdict"] = "PASS" if hits else "NO_QUALIFYING_METRIC"
    observed_missing_contrasts = sorted(
        (item["irregular"], item["control"])
        for item in contrast_status
        if not item["evaluable"]
    )
    if observed_missing_contrasts != sorted(AP37_MISSING_CONTRASTS):
        raise RuntimeError(
            f"A-P3-7 missing-contrast mismatch: {observed_missing_contrasts}"
        )
    h1b_hits = [item for item in waste_contrasts if item["passes_h1b"]]

    controller = load_json(root / "controller-state.json")
    trace_count = sum(
        len(list((Path(summary["run_dir"]) / "traces").glob("*.pt.trace.json*")))
        for summary in summaries.values()
    )
    other_gpu_processes = []
    for stage in complete_stages:
        path = stage / "other-gpu-processes.json"
        if path.exists():
            other_gpu_processes.extend(load_json(path))
    clean_failures = sum(int(item["clean"]["failed"]) for item in summaries.values())
    moderate_rate_ok = []
    for run_id, summary in summaries.items():
        if summary["load"] == "moderate":
            offered = float(summary["clean"]["offered_rps"])
            requested = float(load_json(Path(summary["run_dir"]) / "client/result.json")["request_rate"])
            moderate_rate_ok.append(abs(offered / requested - 1) <= 0.05)
    ratios = []
    for summary in summaries.values():
        ratios.extend(
            [
                summary["waste"]["padding_fraction"],
                summary["waste"]["graph_miss_rate"],
                summary["waste"]["overflow_rate"],
                summary["waste"]["R64"],
                summary["layer1"]["kv_usage_mean"],
                summary["layer1"]["kv_usage_max"],
            ]
        )
    p10_ap36 = ap36_warmup_stability(
        root / "primary/P10-C00-TP2/saturation"
    )
    same_wave_warmup = {}
    for cell in ("P11-C00", "P03-C11"):
        result = load_json(root / f"primary/{cell}/saturation/client/result.json")
        requests = jsonl(root / f"primary/{cell}/saturation/client/requests.jsonl")
        same_wave_warmup[cell] = sum(
            bool(item["success"])
            and 0
            <= float(item["completed_s"])
            < float(result["warmup_seconds"])
            for item in requests
        )
    boundary_marker = load_json(
        Path(summaries["P01-C01-moderate"]["run_dir"]) / "run-complete.json"
    )
    operational_findings = {
        "p10_tp2_non_stabilization": {
            **p10_ap36,
            "status": "PATTERN_CONDITIONED_OPERATIONAL_FINDING",
            "accepted_measurement": False,
            "comparison": (
                "all synthetic pattern runs passed their applicable registered "
                "warm-up gates; orchestrator adjudication"
            ),
            "same_wave_synthetic_warmup_completions": same_wave_warmup,
        },
        "long_context_drain": {
            "run_id": "P10-C01-saturation",
            "drain_seconds": summaries["P10-C01-saturation"]["drain_seconds"],
            "amended_budget_seconds": 600,
            "quarantined": False,
        },
        "failure_boundary": {
            "run_id": "P01-C01-moderate",
            "clean_failures": summaries["P01-C01-moderate"]["clean"]["failed"],
            "excluded_window_failures": boundary_marker["client"][
                "excluded_window_failures"
            ],
            "excluded_failure_kinds": boundary_marker["client"][
                "excluded_window_failure_kinds"
            ],
        },
        "layer2_sampling": {
            "completed_c00_moderate_patterns": len(completed_c00_patterns),
            "evaluable_c00_moderate_patterns": sum(
                operator_share_table[pattern]["moderate"]["status"]
                == "EVALUABLE"
                for pattern in completed_c00_patterns
            ),
            "invalid_windows": len(invalid_primary_windows),
            "classifiable_fraction_min": min(
                window["classifiable_fraction"]
                for windows in layer2.values()
                for window in windows
            ),
        },
    }
    sanity = {
        "numeric": {
            "completed_throughput_rps": numeric_sanity(
                [item["clean"]["completed_throughput_rps"] for item in summaries.values()]
            ),
            "token_efficiency_per_ms": numeric_sanity(
                [item["layer1"]["token_efficiency_per_ms"] for item in summaries.values()]
            ),
            "drain_seconds": numeric_sanity([item["drain_seconds"] for item in summaries.values()]),
            "layer1_clean_steps": numeric_sanity(
                [item["layer1"]["records_clean"] for item in summaries.values()]
            ),
            "operator_classifiable_fraction": numeric_sanity(
                [window["classifiable_fraction"] for windows in layer2.values() for window in windows]
            ),
            "waste_ratios": numeric_sanity(ratios),
            "kendall_tau_b": numeric_sanity([item["tau_b"] for item in tau]),
            "operator_share": numeric_sanity(
                [
                    share
                    for windows in layer2.values()
                    for window in windows
                    for share in window["shares"].values()
                ]
            ),
            "waste_contrast_effect": numeric_sanity(
                [item["point"] for item in waste_contrasts]
            ),
        },
        "invariants": {
            "ap37_primary_runs_40": len(markers) == 40,
            "ap37_confirmation_runs_0": len(confirmation_markers) == 0,
            "accepted_run_ids_unique": len(summaries) == len(markers),
            "complete_cells_20": len(complete_cells) == 20,
            "missing_cells_exact": missing_cells == sorted(AP37_MISSING_CELLS),
            "completed_c00_patterns_9": len(completed_c00_patterns) == 9,
            "clean_duration_240": all(item["clean"]["duration_s"] == 240 for item in summaries.values()),
            "clean_failures_zero": clean_failures == 0,
            "moderate_rate_within_5pct": all(moderate_rate_ok),
            "layer1_footer_invariants": all(
                all(load_json(Path(item["run_dir"]) / "run-complete.json")["layer1"]["invariants"].values())
                for item in summaries.values()
            ),
            "ratios_in_unit_interval": all(0 <= value <= 1 for value in ratios),
            "trace_count_accepted_72": trace_count == 72,
            "missing_trace_count_accepted_8": controller.get("missing_trace_files") == 8,
            "controller_frozen_at_ap36_failure": controller.get("status") == "failed"
            and controller.get("completed_measured_runs") == 40,
            "complete_stage_count_12": len(complete_stages) == 12,
            "clock_and_load_snapshots_complete": all(
                (stage / "clocks-before.txt").exists()
                and (stage / "clocks-after.txt").exists()
                and (stage / "loadavg-before.txt").exists()
                and (stage / "loadavg-after.txt").exists()
                for stage in complete_stages
            ),
            "other_gpu_processes_absent": not other_gpu_processes,
            "drain_quarantine_under_20pct": controller.get("drain_quarantined_runs", 0) / 40 <= 0.20,
            "no_layer2_parser_issues": not layer2_issues,
            "h1b_evaluable_contrasts_6": sum(
                item["evaluable"] for item in contrast_status
            )
            == 6,
            "h1b_missing_contrasts_exact": observed_missing_contrasts
            == sorted(AP37_MISSING_CONTRASTS),
            "ap36_operational_finding_reproduced": math.isclose(
                p10_ap36["normalized_drift"], 0.367639109533929
            )
            and p10_ap36["warmup_completions"] == 17,
            "patterns_not_all_identical_throughput": len(
                {round(item["clean"]["completed_throughput_rps"], 9) for item in summaries.values()}
            )
            > 1,
        },
        "declared_deviations": {
            "missing_saturation_traces": 8,
            "missing_cells": missing_cells,
            "missing_confirmations": ["P10", "P06", "P03", "P01"],
            "unaccepted_canonical_run_markers_excluded": unaccepted_markers,
            "invalid_primary_layer2_windows": invalid_primary_windows,
            "moe_layer_cv": "N/A: layer scopes do not cover >=80% of MoE GEMM time",
        },
    }
    if not all(sanity["invariants"].values()):
        raise RuntimeError(f"data sanity red flag: {sanity['invariants']}")
    return {
        "schema": 1,
        "analysis_seed": SEED,
        "bootstrap_resamples": BOOTSTRAPS,
        "matrix": {
            "primary_runs": len(markers),
            "confirmation_runs": len(confirmation_markers),
            "complete_cells": complete_cells,
            "missing_cells": missing_cells,
            "completed_c00_patterns": completed_c00_patterns,
            "trace_files": trace_count,
            "drain_quarantined_runs": controller.get("drain_quarantined_runs", 0),
            "clean_window_failures": clean_failures,
            "gpu_hours_total": controller["gpu_hours_total"],
        },
        "runs": summaries,
        "all_operator_windows": all_operator_windows,
        "operator_windows": operator_table,
        "operator_share_table": operator_share_table,
        "operator_mode_segments": operator_mode_segments,
        "ranking": {
            "tests": total_ranking_tests,
            "completed_patterns": completed_c00_patterns,
            "evaluable_patterns": [
                pattern
                for pattern in completed_c00_patterns
                if operator_share_table[pattern]["moderate"]["status"]
                == "EVALUABLE"
            ],
            "missing_patterns": ["P05", "P11"],
            "candidates": inversions,
            "accepted_inversions": accepted_inversions,
            "kendall_tau_b": tau,
            "invalid_primary_windows": invalid_primary_windows,
        },
        "waste_contrasts": waste_contrasts,
        "waste_contrast_status": contrast_status,
        "waste_thresholds": thresholds,
        "robustness": {
            "mixed_interference": "leave-one-pattern-out fits applied within config/load",
            "operator_ranking": "two fixed wall-separated windows; no separate LOAO procedure was preregistered",
            "confirmation_runs": "NOT_EVALUABLE: all four confirmations are missing",
        },
        "operational_findings": operational_findings,
        "hypothesis": {
            "H1a": partial_verdict(bool(accepted_inversions)),
            "H1b": partial_verdict(bool(h1b_hits)),
            "compound": "CONFIRMED" if accepted_inversions and h1b_hits else "PARTIAL" if accepted_inversions or h1b_hits else "INCONCLUSIVE",
            "h1b_hits": h1b_hits,
            "refutation_allowed": False,
            "logical_asymmetry": "A-P3-7 permits existential confirmation but incomplete coverage forbids refutation",
        },
        "sanity": sanity,
    }


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--root", type=Path, required=True)
    parser.add_argument("--private", type=Path, required=True)
    parser.add_argument("--out", type=Path, required=True)
    args = parser.parse_args()
    result = analyze(args.root, args.private)
    args.out.parent.mkdir(parents=True, exist_ok=True)
    temporary = args.out.with_suffix(args.out.suffix + ".tmp")
    temporary.write_text(json.dumps(result, indent=2, sort_keys=True, allow_nan=False) + "\n")
    temporary.replace(args.out)
    print(json.dumps({"out": str(args.out), "hypothesis": result["hypothesis"], "sanity": result["sanity"]}, sort_keys=True))


if __name__ == "__main__":
    main()
