scripts/analyze_study5.py

#!/usr/bin/env python3
"""Fail-closed analysis for Study 5 repository-scale harness experiments.

The task is the independent sampling unit.  All contrasts first average over
models within task, then use a task-cluster bootstrap and a Monte-Carlo exact
sign-flip reference distribution.  Raw cells are accepted only when they match
the frozen cell manifests and final experiment reports exactly.
"""

from __future__ import annotations

import argparse
from collections import Counter, defaultdict
import csv
from hashlib import sha256
import json
import math
from pathlib import Path
import random
import statistics
from typing import Any, Callable, Iterable, Sequence

import matplotlib.pyplot as plt
import numpy as np

from agent_harness.specs import load_harnesses, load_tasks


EXPERIMENTS = {"E13": 1440, "E14": 540, "E15": 540, "E16": 306}
REVISIONS = {
    "E13": "7f4de67853deab34aca5a9ceaaf7e9f85b901088",
    "E14": "7f4de67853deab34aca5a9ceaaf7e9f85b901088",
    "E15": "7f4de67853deab34aca5a9ceaaf7e9f85b901088",
    "E16": "6f8c86a1a0737db56c119f90b7e9516344d13dc3",
}
MODELS = ("M002", "M003", "M004")
BOOTSTRAPS = 20_000
BOOTSTRAP_SEED = 20260720
SIGN_FLIPS = 20_000
EPS = 1e-15


class Study5AnalysisError(RuntimeError):
    pass


def canonical_hash(value: Any) -> str:
    return sha256(json.dumps(value, sort_keys=True, separators=(",", ":")).encode()).hexdigest()


def raw_digest(root: Path, experiment: str) -> str:
    digest = sha256()
    for path in sorted((root / "results" / "raw" / experiment).glob("*/*/*/final_metrics.json")):
        digest.update(str(path.relative_to(root)).encode())
        digest.update(b"\0")
        digest.update(path.read_bytes())
        digest.update(b"\0")
    return digest.hexdigest()


def holm(values: Sequence[float]) -> list[float]:
    order = sorted(range(len(values)), key=values.__getitem__)
    adjusted = [1.0] * len(values)
    running = 0.0
    for rank, index in enumerate(order):
        running = max(running, min(1.0, (len(values) - rank) * values[index]))
        adjusted[index] = running
    return adjusted


def cluster_bootstrap(values: Sequence[float], seed_offset: int = 0) -> tuple[float, float]:
    if not values:
        return math.nan, math.nan
    generator = random.Random(BOOTSTRAP_SEED + seed_offset)
    n = len(values)
    draws = sorted(
        statistics.fmean(values[generator.randrange(n)] for _ in range(n))
        for _ in range(BOOTSTRAPS)
    )
    return draws[int(0.025 * BOOTSTRAPS)], draws[int(0.975 * BOOTSTRAPS) - 1]


def sign_flip(values: Sequence[float], seed_offset: int = 0) -> float:
    if not values or all(abs(value) < EPS for value in values):
        return 1.0
    observed = abs(statistics.fmean(values))
    generator = random.Random(BOOTSTRAP_SEED + 100_000 + seed_offset)
    extreme = 0
    for _ in range(SIGN_FLIPS):
        value = abs(statistics.fmean(
            item if generator.getrandbits(1) else -item for item in values
        ))
        extreme += value + EPS >= observed
    return (extreme + 1) / (SIGN_FLIPS + 1)


def exact_mcnemar(left: Sequence[int], right: Sequence[int]) -> tuple[int, int, float]:
    n10 = sum(a == 1 and b == 0 for a, b in zip(left, right))
    n01 = sum(a == 0 and b == 1 for a, b in zip(left, right))
    n = n10 + n01
    if not n:
        return n10, n01, 1.0
    tail = sum(math.comb(n, k) for k in range(min(n10, n01) + 1)) / 2**n
    return n10, n01, min(1.0, 2 * tail)


def contrast_record(name: str, task_effects: Sequence[float], seed_offset: int = 0) -> dict[str, Any]:
    low, high = cluster_bootstrap(task_effects, seed_offset)
    return {
        "contrast": name,
        "independent_tasks": len(task_effects),
        "estimate": statistics.fmean(task_effects),
        "cluster_bootstrap_ci_low": low,
        "cluster_bootstrap_ci_high": high,
        "sign_flip_p": sign_flip(task_effects, seed_offset),
    }


def _mechanism(trajectory: Path, gold_files: set[str]) -> dict[str, int]:
    retrieved: set[str] = set()
    read: set[str] = set()
    accepted = False
    for line in trajectory.read_text(encoding="utf-8").splitlines():
        event = json.loads(line)
        payload = event.get("payload", {})
        if event.get("event_type") == "edit" and payload.get("accepted"):
            accepted = True
        if accepted:
            continue
        if event.get("event_type") == "retrieval_candidate" and payload.get("path"):
            retrieved.add(str(payload["path"]))
        if event.get("event_type") == "file_read" and payload.get("path"):
            read.add(str(payload["path"]))
    return {
        "gold_retrieved_before_edit": int(bool(retrieved & gold_files)),
        "gold_read_before_edit": int(bool(read & gold_files)),
        "unique_retrieved_before_edit": len(retrieved),
        "unique_read_before_edit": len(read),
    }


def _complete_report(root: Path, experiment: str, expected: int) -> tuple[Path, dict[str, Any]]:
    candidates: list[tuple[Path, dict[str, Any]]] = []
    for path in sorted((root / "results" / "reports").glob(f"{experiment}_*.json")):
        if path.name.endswith("_progress.json"):
            continue
        report = json.loads(path.read_text(encoding="utf-8"))
        if report.get("run_count") == expected:
            candidates.append((path, report))
    if len(candidates) != 1:
        raise Study5AnalysisError(
            f"{experiment} requires exactly one complete final report; found {len(candidates)}"
        )
    return candidates[0]


def discover(root: Path) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    tasks = load_tasks(root)
    all_rows: list[dict[str, Any]] = []
    audit: dict[str, Any] = {"schema_version": 1, "experiments": {}}
    selection = json.loads((root / "configs" / "study5" / "E16_selection.json").read_text())
    selection_check = dict(selection)
    expected_selection_hash = selection_check.pop("selection_sha256")
    if canonical_hash(selection_check) != expected_selection_hash:
        raise Study5AnalysisError("E16 selection hash mismatch")

    for experiment, expected_count in EXPERIMENTS.items():
        manifest_path = root / "configs" / "study5" / f"{experiment}_cells.json"
        frozen = json.loads(manifest_path.read_text(encoding="utf-8"))
        frozen_check = dict(frozen)
        expected_design_hash = frozen_check.pop("design_sha256")
        if canonical_hash(frozen_check) != expected_design_hash:
            raise Study5AnalysisError(f"{experiment} frozen design hash mismatch")
        expected = {
            (cell["task_id"], cell["harness_id"], cell["interface_id"], cell["model_id"])
            for cell in frozen["cells"]
        }
        if len(expected) != expected_count or frozen["planned_cells"] != expected_count:
            raise Study5AnalysisError(f"{experiment} frozen grid count mismatch")
        final_paths = sorted((root / "results" / "raw" / experiment).glob("*/*/*/final_metrics.json"))
        if len(final_paths) != expected_count:
            raise Study5AnalysisError(
                f"{experiment} requires {expected_count} final metrics, observed {len(final_paths)}"
            )
        seen: set[tuple[str, str, str, str]] = set()
        revisions: Counter[str] = Counter()
        for final_path in final_paths:
            directory = final_path.parent
            required = [
                directory / "run_manifest.json", directory / "trajectory.jsonl",
                directory / "messages.json", directory / "model.patch", directory / "validation.json",
            ]
            if not all(path.is_file() for path in required):
                raise Study5AnalysisError(f"incomplete raw cell: {directory}")
            final = json.loads(final_path.read_text(encoding="utf-8"))
            run_manifest = json.loads(required[0].read_text(encoding="utf-8"))
            identity = run_manifest["identity"]
            key = (
                identity["task_id"], final["retrieval_harness_id"],
                final["edit_interface_id"], identity["model_id"],
            )
            if key not in expected or key in seen:
                raise Study5AnalysisError(f"unexpected or duplicate {experiment} identity: {key}")
            if identity["harness_id"] != f"{key[1]}__{key[2]}":
                raise Study5AnalysisError(f"composite treatment mismatch: {directory}")
            if identity["context_budget"] != 65536 or identity["seed"] != 0:
                raise Study5AnalysisError(f"inference identity drift: {directory}")
            if any(len(item.get("after_instances", [])) != 1 for item in final["residency_transitions"]):
                raise Study5AnalysisError(f"non-exclusive model residency: {directory}")
            seen.add(key)
            revisions[identity["code_revision"]] += 1
            task = tasks[identity["task_id"]]
            mechanism = _mechanism(required[1], set(task.gold_files))
            all_rows.append({
                "experiment_id": experiment,
                "run_id": run_manifest["run_id"],
                "task_id": identity["task_id"],
                "repository_sha": identity["repository_sha"],
                "model_id": identity["model_id"],
                "harness_id": final["retrieval_harness_id"],
                "interface_id": final["edit_interface_id"],
                "resolved_at_1": int(bool(final["resolved_at_1"])),
                "accepted_edit_cell": int(bool(final["accepted_edit_cell"])),
                "applicable_final_patch": int(bool(final["applicable_final_patch"])),
                "exact_modified_file_match": int(bool(final["exact_modified_file_match"])),
                "fail_to_pass": int(bool(final["fail_to_pass"])),
                "search_gold_any": int(final["search_localization_metrics"]["file_recall_at_10"] > 0),
                "read_gold_any": int(final["read_localization_metrics"]["file_recall_at_10"] > 0),
                "search_mrr": float(final["search_localization_metrics"]["mrr"]),
                "read_mrr": float(final["read_localization_metrics"]["mrr"]),
                **mechanism,
                "total_tokens": int(final["usage"]["total_tokens"]),
                "prompt_tokens": int(final["usage"]["prompt_tokens"]),
                "completion_tokens": int(final["usage"]["completion_tokens"]),
                "elapsed_seconds": float(final["elapsed_seconds"]),
                "model_elapsed_seconds": float(final["model_elapsed_seconds"]),
                "model_switch_seconds": float(final["model_switch_seconds"]),
                "model_switch_count": int(final["model_switch_count"]),
                "tool_calls": int(final["tool_calls"]),
                "protocol_violation_count": len(final["protocol_violations"]),
                "failure_stage": str(final["failure_stage"]),
                "finished_reason": str(final["finished_reason"]),
            })
        if seen != expected or set(revisions) != {REVISIONS[experiment]}:
            raise Study5AnalysisError(f"{experiment} grid/revision mismatch: {revisions}")
        report_path, report = _complete_report(root, experiment, expected_count)
        counts = {
            "run_count": expected_count,
            "accepted_edit_count": sum(row["accepted_edit_cell"] for row in all_rows if row["experiment_id"] == experiment),
            "applicable_patch_count": sum(row["applicable_final_patch"] for row in all_rows if row["experiment_id"] == experiment),
            "resolved_count": sum(row["resolved_at_1"] for row in all_rows if row["experiment_id"] == experiment),
        }
        if report.get("code_revision") != REVISIONS[experiment]:
            raise Study5AnalysisError(f"{experiment} final report revision mismatch")
        if any(report.get(key) != value for key, value in counts.items()):
            raise Study5AnalysisError(f"{experiment} raw ledger/final report count mismatch")
        if report.get("manifest_sha256") != expected_design_hash:
            raise Study5AnalysisError(f"{experiment} report/design manifest hash mismatch")
        digest = raw_digest(root, experiment)
        if experiment in ("E13", "E15") and digest != selection["source_raw_sha256"][experiment]:
            raise Study5AnalysisError(f"{experiment} raw hash no longer matches E16 screening freeze")
        audit["experiments"][experiment] = {
            "cells": expected_count,
            "unique_identities": len(seen),
            "code_revision": REVISIONS[experiment],
            "design_sha256": expected_design_hash,
            "raw_final_metrics_sha256": digest,
            "final_report": str(report_path.relative_to(root)),
            "final_report_sha256": sha256(report_path.read_bytes()).hexdigest(),
            **counts,
        }
    audit["total_cells"] = len(all_rows)
    audit["selection_sha256"] = expected_selection_hash
    audit["passed"] = len(all_rows) == sum(EXPERIMENTS.values())
    audit["audit_sha256"] = canonical_hash(audit)
    return all_rows, audit


def summarize(rows: Sequence[dict[str, Any]], keys: Sequence[str]) -> list[dict[str, Any]]:
    grouped: dict[tuple[Any, ...], list[dict[str, Any]]] = defaultdict(list)
    for row in rows:
        grouped[tuple(row[key] for key in keys)].append(row)
    output = []
    for group, values in sorted(grouped.items()):
        record = dict(zip(keys, group))
        record.update({
            "cells": len(values),
            "tasks": len({row["task_id"] for row in values}),
            "search_gold_rate": statistics.fmean(row["search_gold_any"] for row in values),
            "read_gold_rate": statistics.fmean(row["read_gold_any"] for row in values),
            "accepted_edit_rate": statistics.fmean(row["accepted_edit_cell"] for row in values),
            "applicable_patch_rate": statistics.fmean(row["applicable_final_patch"] for row in values),
            "resolved_rate": statistics.fmean(row["resolved_at_1"] for row in values),
            "mean_total_tokens": statistics.fmean(row["total_tokens"] for row in values),
            "mean_wall_seconds": statistics.fmean(row["elapsed_seconds"] for row in values),
            "mean_model_switch_seconds": statistics.fmean(row["model_switch_seconds"] for row in values),
            "mean_protocol_violations": statistics.fmean(row["protocol_violation_count"] for row in values),
        })
        output.append(record)
    return output


def _task_factor_effects(
    rows: Sequence[dict[str, Any]], endpoint: str, contrast: Callable[[dict[str, int]], int]
) -> list[float]:
    grouped: dict[str, dict[int, list[float]]] = defaultdict(lambda: defaultdict(list))
    for row in rows:
        code = {"L": int(row["L"]), "S": int(row["S"]), "D": int(row["D"])}
        grouped[row["task_id"]][contrast(code)].append(float(row[endpoint]))
    effects = []
    for task in sorted(grouped):
        if set(grouped[task]) != {-1, 1}:
            raise Study5AnalysisError(f"incomplete factorial contrast for task {task}")
        effects.append(statistics.fmean(grouped[task][1]) - statistics.fmean(grouped[task][-1]))
    return effects


def e13_factorial(rows: Sequence[dict[str, Any]], harness_specs: dict[str, Any]) -> dict[str, Any]:
    selected = [dict(row) for row in rows if row["experiment_id"] == "E13"]
    for row in selected:
        spec = harness_specs[row["harness_id"]]
        row.update(L=int(spec.lexical), S=int(spec.syntax == "tree_sitter"), D=int(spec.dense))
    contrasts: dict[str, Callable[[dict[str, int]], int]] = {
        "L": lambda x: 1 if x["L"] else -1,
        "S": lambda x: 1 if x["S"] else -1,
        "D": lambda x: 1 if x["D"] else -1,
        "LxS": lambda x: (1 if x["L"] else -1) * (1 if x["S"] else -1),
        "LxD": lambda x: (1 if x["L"] else -1) * (1 if x["D"] else -1),
        "SxD": lambda x: (1 if x["S"] else -1) * (1 if x["D"] else -1),
        "LxSxD": lambda x: (1 if x["L"] else -1) * (1 if x["S"] else -1) * (1 if x["D"] else -1),
    }
    endpoints = ("resolved_at_1", "accepted_edit_cell", "applicable_final_patch", "search_gold_any", "read_gold_any")
    output: dict[str, Any] = {"estimand": "marginal risk difference averaged over models and remaining factors", "endpoints": {}}
    offset = 0
    for endpoint in endpoints:
        records = []
        for name, function in contrasts.items():
            effects = _task_factor_effects(selected, endpoint, function)
            records.append(contrast_record(f"E13_{name}_{endpoint}", effects, offset))
            offset += 1
        adjusted = holm([record["sign_flip_p"] for record in records])
        for record, value in zip(records, adjusted):
            record["sign_flip_p_holm_within_endpoint"] = value
        output["endpoints"][endpoint] = records
    model_interactions = []
    for endpoint in ("resolved_at_1", "accepted_edit_cell", "applicable_final_patch"):
        for factor in ("L", "S", "D"):
            function = contrasts[factor]
            per_model: dict[str, list[float]] = {}
            for model in MODELS:
                per_model[model] = _task_factor_effects(
                    [row for row in selected if row["model_id"] == model], endpoint, function
                )
            for left, right in (("M002", "M003"), ("M002", "M004"), ("M003", "M004")):
                effects = [a - b for a, b in zip(per_model[left], per_model[right])]
                model_interactions.append(contrast_record(
                    f"E13_{factor}_{endpoint}_{left}_minus_{right}", effects, offset
                ))
                offset += 1
    adjusted = holm([record["sign_flip_p"] for record in model_interactions])
    for record, value in zip(model_interactions, adjusted):
        record["sign_flip_p_holm"] = value
    output["model_interactions"] = model_interactions
    return output


def _cell_lookup(rows: Sequence[dict[str, Any]], dimensions: Sequence[str]) -> dict[tuple[Any, ...], dict[str, Any]]:
    lookup = {tuple(row[key] for key in dimensions): row for row in rows}
    if len(lookup) != len(rows):
        raise Study5AnalysisError(f"duplicate cells for lookup dimensions {dimensions}")
    return lookup


def e14_interactions(rows: Sequence[dict[str, Any]]) -> dict[str, Any]:
    selected = [row for row in rows if row["experiment_id"] == "E14"]
    lookup = _cell_lookup(selected, ("task_id", "model_id", "harness_id", "interface_id"))
    tasks = sorted({row["task_id"] for row in selected})
    retrieval_pairs = (("H006", "H000"), ("H007", "H000"), ("H007", "H006"))
    action_pairs = (("P002", "P001"), ("P003", "P001"), ("P003", "P002"))
    output: dict[str, Any] = {"estimand": "task-level difference-in-differences averaged over models", "endpoints": {}}
    offset = 100
    for endpoint in ("accepted_edit_cell", "applicable_final_patch", "resolved_at_1"):
        records = []
        for high_h, low_h in retrieval_pairs:
            for high_p, low_p in action_pairs:
                effects = []
                for task in tasks:
                    model_effects = []
                    for model in MODELS:
                        def value(h: str, p: str) -> float:
                            return float(lookup[(task, model, h, p)][endpoint])
                        model_effects.append((value(high_h, high_p) - value(low_h, high_p)) -
                                             (value(high_h, low_p) - value(low_h, low_p)))
                    effects.append(statistics.fmean(model_effects))
                records.append(contrast_record(
                    f"E14_{high_h}-{low_h}_x_{high_p}-{low_p}_{endpoint}", effects, offset
                ))
                offset += 1
        adjusted = holm([record["sign_flip_p"] for record in records])
        for record, value in zip(records, adjusted):
            record["sign_flip_p_holm_within_endpoint"] = value
        output["endpoints"][endpoint] = records
    return output


def _paired_harness_effects(
    selected: Sequence[dict[str, Any]], left: str, right: str, endpoint: str
) -> list[float]:
    lookup = _cell_lookup(selected, ("task_id", "model_id", "harness_id"))
    tasks = sorted({row["task_id"] for row in selected})
    return [statistics.fmean(
        float(lookup[(task, model, left)][endpoint]) - float(lookup[(task, model, right)][endpoint])
        for model in MODELS
    ) for task in tasks]


def e15_blocks(rows: Sequence[dict[str, Any]]) -> dict[str, Any]:
    selected = [row for row in rows if row["experiment_id"] == "E15"]
    blocks = {
        "graph": (("H007", "H008"), ("H009", "H008")),
        "query_policy": (("H010", "H008"),),
        "search_interface": (("H011", "H008"),),
        "query_by_interface": (("H010", "H008"), ("H011", "H008"), ("H012", "H008")),
        "packing": (("H013", "H008"), ("H014", "H008"), ("H015", "H008")),
    }
    endpoints = (
        "resolved_at_1", "accepted_edit_cell", "applicable_final_patch",
        "gold_retrieved_before_edit", "gold_read_before_edit", "total_tokens",
        "elapsed_seconds", "model_switch_seconds", "protocol_violation_count",
    )
    output: dict[str, Any] = {}
    offset = 200
    for block, pairs in blocks.items():
        block_records = []
        for endpoint in endpoints:
            endpoint_records = []
            for left, right in pairs:
                effects = _paired_harness_effects(selected, left, right, endpoint)
                endpoint_records.append(contrast_record(
                    f"E15_{block}_{left}_minus_{right}_{endpoint}", effects, offset
                ))
                offset += 1
            adjusted = holm([record["sign_flip_p"] for record in endpoint_records])
            for record, value in zip(endpoint_records, adjusted):
                record["sign_flip_p_holm_within_block_endpoint"] = value
            block_records.extend(endpoint_records)
        output[block] = block_records
    return output


def rankdata(values: dict[str, float], higher_better: bool = True) -> dict[str, float]:
    ordered = sorted(values, key=lambda key: (-values[key] if higher_better else values[key], key))
    output: dict[str, float] = {}
    index = 0
    while index < len(ordered):
        end = index + 1
        while end < len(ordered) and abs(values[ordered[end]] - values[ordered[index]]) < EPS:
            end += 1
        average_rank = (index + 1 + end) / 2
        for key in ordered[index:end]:
            output[key] = average_rank
        index = end
    return output


def spearman(left: dict[str, float], right: dict[str, float], higher_better: bool = True) -> float:
    keys = sorted(set(left) & set(right))
    a = rankdata({key: left[key] for key in keys}, higher_better)
    b = rankdata({key: right[key] for key in keys}, higher_better)
    av = statistics.fmean(a.values()); bv = statistics.fmean(b.values())
    numerator = sum((a[key] - av) * (b[key] - bv) for key in keys)
    denominator = math.sqrt(sum((a[key] - av) ** 2 for key in keys) * sum((b[key] - bv) ** 2 for key in keys))
    return numerator / denominator if denominator else math.nan


def pareto(summary: Sequence[dict[str, Any]]) -> list[str]:
    def dominates(a: dict[str, Any], b: dict[str, Any]) -> bool:
        av = (a["resolved_rate"], a["applicable_patch_rate"], a["accepted_edit_rate"],
              -a["mean_total_tokens"], -a["mean_wall_seconds"])
        bv = (b["resolved_rate"], b["applicable_patch_rate"], b["accepted_edit_rate"],
              -b["mean_total_tokens"], -b["mean_wall_seconds"])
        return all(x >= y for x, y in zip(av, bv)) and any(x > y for x, y in zip(av, bv))
    return sorted(row["harness_id"] for row in summary if not any(
        other["harness_id"] != row["harness_id"] and dominates(other, row) for other in summary
    ))


def e16_validation(rows: Sequence[dict[str, Any]], selection: dict[str, Any]) -> dict[str, Any]:
    selected = [row for row in rows if row["experiment_id"] == "E16"]
    summary = summarize(selected, ("harness_id",))
    by_harness = {row["harness_id"]: row for row in summary}
    screening = selection["metrics"]
    metric_pairs = {
        "resolved": ("resolved_rate", "resolved_rate", True),
        "accepted": ("accepted_edit_rate", "accepted_edit_rate", True),
        "applicable": ("applicable_patch_rate", "applicable_patch_rate", True),
        "tokens": ("mean_total_tokens", "mean_total_tokens", False),
        "wall_seconds": ("mean_wall_seconds", "mean_wall_seconds", False),
    }
    rank_stability = {}
    for name, (screen_key, held_key, high) in metric_pairs.items():
        left = {key: float(screening[key][screen_key]) for key in by_harness}
        right = {key: float(by_harness[key][held_key]) for key in by_harness}
        rank_stability[name] = {
            "screening": left, "held_out": right,
            "spearman_rank_correlation": spearman(left, right, high),
        }
    paired = []
    offset = 400
    for harness in sorted(by_harness):
        if harness == "H000":
            continue
        for endpoint in ("resolved_at_1", "accepted_edit_cell", "applicable_final_patch", "total_tokens", "elapsed_seconds"):
            effects = _paired_harness_effects(selected, harness, "H000", endpoint)
            paired.append(contrast_record(f"E16_{harness}_minus_H000_{endpoint}", effects, offset))
            offset += 1
    adjusted = holm([row["sign_flip_p"] for row in paired if row["contrast"].endswith("accepted_edit_cell")])
    for row, value in zip([r for r in paired if r["contrast"].endswith("accepted_edit_cell")], adjusted):
        row["sign_flip_p_holm_accepted_family"] = value
    return {
        "harness_summary": summary,
        "rank_stability": rank_stability,
        "paired_vs_exact_baseline": paired,
        "held_out_pareto_frontier": pareto(summary),
        "resolution_event_warning": "Only one E16 cell resolved; resolution ranks are descriptive and not evidence of equivalence.",
    }


def failure_taxonomy(rows: Sequence[dict[str, Any]]) -> list[dict[str, Any]]:
    grouped: Counter[tuple[str, str, str]] = Counter(
        (row["experiment_id"], row["model_id"], row["failure_stage"]) for row in rows
    )
    return [
        {"experiment_id": experiment, "model_id": model, "failure_stage": stage, "cells": count}
        for (experiment, model, stage), count in sorted(grouped.items())
    ]


def write_csv(path: Path, rows: Sequence[dict[str, Any]]) -> None:
    if not rows:
        return
    keys = list(rows[0])
    with path.open("w", encoding="utf-8", newline="") as handle:
        writer = csv.DictWriter(
            handle, fieldnames=keys, extrasaction="ignore", lineterminator="\n"
        )
        writer.writeheader(); writer.writerows(rows)


def figures(output: Path, summaries: Sequence[dict[str, Any]], e16: dict[str, Any]) -> list[str]:
    files: list[str] = []
    e13 = [row for row in summaries if row["experiment_id"] == "E13"]
    harnesses = [row["harness_id"] for row in e13]
    x = np.arange(len(harnesses))
    fig, ax = plt.subplots(figsize=(10, 5.4))
    ax.plot(x, [row["search_gold_rate"] for row in e13], "o-", label="Gold retrieved")
    ax.plot(x, [row["accepted_edit_rate"] for row in e13], "s-", label="Accepted edit")
    ax.plot(x, [row["applicable_patch_rate"] for row in e13], "^-", label="Applicable patch")
    ax.plot(x, [row["resolved_rate"] for row in e13], "D-", label="Resolved@1")
    ax.set_xticks(x, harnesses, rotation=45); ax.set_ylim(-.02, 1.02)
    ax.set_ylabel("Cell rate"); ax.set_title("E13 stage-aware outcome funnel by retrieval harness")
    ax.grid(axis="y", alpha=.25); ax.legend(ncol=2); fig.tight_layout()
    for suffix in ("png", "pdf"):
        path = output / f"figure_e13_funnel.{suffix}"
        metadata = {"CreationDate": None, "ModDate": None} if suffix == "pdf" else None
        fig.savefig(path, dpi=220, metadata=metadata); files.append(path.name)
    plt.close(fig)

    held = e16["harness_summary"]
    fig, ax = plt.subplots(figsize=(7.5, 5.5))
    for row in held:
        ax.scatter(row["mean_total_tokens"], row["applicable_patch_rate"], s=70)
        ax.annotate(row["harness_id"], (row["mean_total_tokens"], row["applicable_patch_rate"]), xytext=(4, 4), textcoords="offset points")
    ax.set_xlabel("Mean total tokens (lower is better)"); ax.set_ylabel("Applicable-patch rate")
    ax.set_title("E16 held-out quality–cost behavior"); ax.grid(alpha=.25); fig.tight_layout()
    for suffix in ("png", "pdf"):
        path = output / f"figure_e16_quality_cost.{suffix}"
        metadata = {"CreationDate": None, "ModDate": None} if suffix == "pdf" else None
        fig.savefig(path, dpi=220, metadata=metadata); files.append(path.name)
    plt.close(fig)
    return files


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--root", type=Path, default=Path(__file__).resolve().parents[1])
    args = parser.parse_args()
    root = args.root.resolve()
    rows, audit = discover(root)
    harness_specs = load_harnesses(root)
    selection = json.loads((root / "configs" / "study5" / "E16_selection.json").read_text())
    summaries = summarize(rows, ("experiment_id", "harness_id"))
    model_summaries = summarize(rows, ("experiment_id", "model_id"))
    analysis = {
        "schema_version": 1,
        "bootstrap_draws": BOOTSTRAPS,
        "sign_flip_draws": SIGN_FLIPS,
        "seed": BOOTSTRAP_SEED,
        "independent_unit": "task",
        "e13_factorial": e13_factorial(rows, harness_specs),
        "e14_retrieval_action": e14_interactions(rows),
        "e15_navigation_packing": e15_blocks(rows),
        "e16_held_out": e16_validation(rows, selection),
        "failure_taxonomy": failure_taxonomy(rows),
        "claim_boundary": "Zero or sparse resolution differences are no evidence of improvement, not evidence of equivalence.",
    }
    output = root / "results" / "derived" / "study5"
    output.mkdir(parents=True, exist_ok=True)
    write_csv(output / "all_cells.csv", rows)
    write_csv(output / "harness_summary.csv", summaries)
    write_csv(output / "model_summary.csv", model_summaries)
    write_csv(output / "failure_taxonomy.csv", analysis["failure_taxonomy"])
    (output / "artifact_audit.json").write_text(json.dumps(audit, indent=2, sort_keys=True) + "\n")
    (output / "statistical_results.json").write_text(json.dumps(analysis, indent=2, sort_keys=True, allow_nan=False) + "\n")
    figure_files = figures(output, summaries, analysis["e16_held_out"])
    manifest = {
        "schema_version": 1,
        "input_audit_sha256": audit["audit_sha256"],
        "analysis_sha256": sha256((output / "statistical_results.json").read_bytes()).hexdigest(),
        "cells_csv_sha256": sha256((output / "all_cells.csv").read_bytes()).hexdigest(),
        "figure_files": figure_files,
        "files": {},
    }
    for path in sorted(output.iterdir()):
        if path.is_file() and path.name != "analysis_manifest.json":
            manifest["files"][path.name] = sha256(path.read_bytes()).hexdigest()
    manifest["manifest_sha256"] = canonical_hash(manifest)
    (output / "analysis_manifest.json").write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n")
    print(json.dumps({
        "passed": audit["passed"], "cells": len(rows), "output": str(output),
        "audit_sha256": audit["audit_sha256"], "analysis_sha256": manifest["analysis_sha256"],
    }, indent=2))


if __name__ == "__main__":
    main()