scripts/analyze_study3.py

#!/usr/bin/env python3
"""Prespecified fail-closed analysis for the frozen E09 protocol study."""

from __future__ import annotations

import argparse
from collections import Counter
import csv
from hashlib import sha256
import json
import math
from pathlib import Path
import random
import statistics
import subprocess
import tomllib
from typing import Any, Sequence

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


MODELS = ("M002", "M003", "M004")
INTERFACES = ("P001", "P002", "P003")
REPOSITORIES = ("R001", "R002", "R003")
BOOTSTRAPS = 20_000
BOOTSTRAP_SEED = 20260719
PRIMARY = ("P1_M003_P002_vs_P001_accepted", "accepted_edit_cell", "M003", "P002", "P001")


class Study3AnalysisError(RuntimeError):
    pass


def percentile(values: Sequence[float], probability: float) -> float:
    ordered = sorted(float(value) for value in values)
    if not ordered:
        raise ValueError("percentile requires values")
    position = (len(ordered) - 1) * probability
    low, high = math.floor(position), math.ceil(position)
    if low == high:
        return ordered[low]
    return ordered[low] * (high - position) + ordered[high] * (position - low)


def bootstrap_mean_ci(
    values: Sequence[float], rng: random.Random, samples: int = BOOTSTRAPS
) -> tuple[float, float]:
    draws = [statistics.fmean(rng.choice(values) for _ in values) for _ in range(samples)]
    return percentile(draws, 0.025), percentile(draws, 0.975)


def exact_mcnemar(left: Sequence[int], right: Sequence[int]) -> tuple[int, int, float]:
    if len(left) != len(right):
        raise ValueError("paired vectors differ in length")
    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))
    discordant = n10 + n01
    if discordant == 0:
        return n10, n01, 1.0
    tail = sum(math.comb(discordant, k) for k in range(min(n10, n01) + 1)) / 2**discordant
    return n10, n01, min(1.0, 2 * tail)


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


def sha256_file(path: Path) -> str:
    value = sha256()
    with path.open("rb") as handle:
        for block in iter(lambda: handle.read(1024 * 1024), b""):
            value.update(block)
    return value.hexdigest()


def git_output(root: Path, *arguments: str) -> str:
    return subprocess.run(
        ["git", *arguments], cwd=root, check=True, capture_output=True, text=True
    ).stdout.strip()


def task_metadata(root: Path) -> dict[str, dict[str, Any]]:
    repositories: dict[str, dict[str, Any]] = {}
    urls: dict[str, str] = {}
    for path in sorted((root / "configs/repositories").glob("R*.toml")):
        value = tomllib.loads(path.read_text(encoding="utf-8"))
        repository_id = str(value["repository_id"])
        repositories[repository_id] = value
        urls[str(value["repository_url"])] = repository_id
    split = [
        line.strip()
        for line in (root / "tasks/splits/study3_protocol.txt").read_text().splitlines()
        if line.strip() and not line.lstrip().startswith("#")
    ]
    result: dict[str, dict[str, Any]] = {}
    for task_id in split:
        value = tomllib.loads(
            (root / "tasks/manifests" / f"{task_id}.toml").read_text(encoding="utf-8")
        )
        repository_id = urls[str(value["repository_url"])]
        result[task_id] = {**value, "repository_id": repository_id}
    if len(result) != 60 or Counter(v["repository_id"] for v in result.values()) != Counter(
        {"R001": 20, "R002": 20, "R003": 20}
    ):
        raise Study3AnalysisError("E09 task metadata is not the frozen 20/20/20 split")
    return result


def _required_artifacts(directory: Path) -> None:
    required = (
        "run_manifest.json",
        "final_metrics.json",
        "messages.json",
        "model.patch",
        "validation.json",
        "trajectory.jsonl",
    )
    missing = [item for item in required if not (directory / item).is_file()]
    if missing:
        raise Study3AnalysisError(f"missing artifacts {missing} in {directory}")


def discover(root: Path) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    tasks = task_metadata(root)
    paths = sorted((root / "results/raw/E09").rglob("final_metrics.json"))
    if len(paths) != 540:
        raise Study3AnalysisError(f"E09 requires 540 final metrics; observed {len(paths)}")
    expected = {
        (task_id, interface, model, 0, 0)
        for task_id in tasks
        for interface in INTERFACES
        for model in MODELS
    }
    identities: set[tuple[Any, ...]] = set()
    run_ids: set[str] = set()
    revisions: Counter[str] = Counter()
    model_keys: Counter[str] = Counter()
    responses = 0
    raw_digest = sha256()
    rows: list[dict[str, Any]] = []
    for path in paths:
        directory = path.parent
        _required_artifacts(directory)
        final = json.loads(path.read_text(encoding="utf-8"))
        manifest = json.loads((directory / "run_manifest.json").read_text(encoding="utf-8"))
        identity = manifest["identity"]
        key = tuple(
            identity[item]
            for item in ("task_id", "harness_id", "model_id", "seed", "repetition")
        )
        if key in identities:
            raise Study3AnalysisError(f"duplicate E09 identity: {key}")
        identities.add(key)
        if key not in expected:
            raise Study3AnalysisError(f"unexpected E09 identity: {key}")
        run_id = str(final["run_id"])
        if run_id in run_ids or run_id != manifest["run_id"]:
            raise Study3AnalysisError(f"duplicate/mismatched run ID: {run_id}")
        run_ids.add(run_id)
        task_id, interface, model_id, seed, repetition = key
        if final["edit_interface_id"] != interface or final["model_id"] != model_id:
            raise Study3AnalysisError(f"manifest/final mismatch in {directory}")
        resolved_model = manifest["resolved_model"]
        runtime = resolved_model["agent_runtime"]
        if runtime["inference_key"] != identity["model_key"]:
            raise Study3AnalysisError(f"inference identity mismatch in {directory}")
        loaded = runtime["native_record"].get("loaded_instances") or []
        if len(loaded) != 1 or int(loaded[0]["config"]["context_length"]) != 65536:
            raise Study3AnalysisError(f"runtime context/residency mismatch in {directory}")
        if any(len(item.get("after_instances", [])) != 1 for item in final["residency_transitions"]):
            raise Study3AnalysisError(f"non-exclusive residency in {directory}")
        response_paths = sorted(directory.glob("model_response_*.json"))
        if len(response_paths) != int(final["model_calls"]):
            raise Study3AnalysisError(f"response count mismatch in {directory}")
        responses += len(response_paths)
        revisions[str(identity["code_revision"])] += 1
        model_keys[str(identity["model_key"])] += 1
        usage = final.get("usage") or {}
        metadata = tasks[str(task_id)]
        row = {
            "run_id": run_id,
            "task_id": str(task_id),
            "repository_id": metadata["repository_id"],
            "language": metadata["language"],
            "difficulty": metadata["difficulty"],
            "interface_id": str(interface),
            "edit_tool": str(final["edit_tool"]),
            "model_id": str(model_id),
            "model_key": str(identity["model_key"]),
            "model_config_hash": str(identity["model_config_hash"]),
            "interface_config_hash": str(identity["harness_hash"]),
            "code_revision": str(identity["code_revision"]),
            "seed": int(seed),
            "repetition": int(repetition),
            "accepted_edit_cell": int(bool(final["accepted_edit_cell"])),
            "edit_attempts": int(final["edit_attempts"]),
            "edit_acceptances": int(final["edit_acceptances"]),
            "applicable_final_patch": int(bool(final["applicable_final_patch"])),
            "resolved_at_1": int(bool(final["resolved_at_1"])),
            "fail_to_pass": int(bool(final["fail_to_pass"])),
            "pass_to_pass": int(bool(final["pass_to_pass"])),
            "exact_modified_file_match": int(bool(final["exact_modified_file_match"])),
            "gold_file_modified_recall": float(final["gold_file_modified_recall"]),
            "read_file_recall_at_10": float(
                final["read_localization_metrics"]["file_recall_at_10"]
            ),
            "protocol_violation_count": len(final["protocol_violations"]),
            "model_calls": int(final["model_calls"]),
            "tool_calls": int(final["tool_calls"]),
            "test_runs": int(final["test_runs"]),
            "prompt_tokens": int(usage.get("prompt_tokens", 0)),
            "completion_tokens": int(usage.get("completion_tokens", 0)),
            "total_tokens": int(usage.get("total_tokens", 0)),
            "elapsed_seconds": float(final["elapsed_seconds"]),
            "finished_reason": str(final["finished_reason"]),
            "failure_stage": str(final["failure_stage"]),
        }
        rows.append(row)
        for artifact in (directory / "run_manifest.json", path):
            raw_digest.update(
                artifact.relative_to(root).as_posix().encode()
                + b"\0"
                + artifact.read_bytes()
                + b"\0"
            )
    if identities != expected:
        raise Study3AnalysisError(f"E09 grid mismatch; missing={sorted(expected-identities)[:5]}")
    if len(revisions) != 1:
        raise Study3AnalysisError(f"mixed E09 revisions: {dict(revisions)}")
    revision = next(iter(revisions))
    preflight_path = root / "results/reports/study3_preflight.json"
    if not preflight_path.is_file():
        raise Study3AnalysisError("Study 3 preflight report is missing")
    preflight = json.loads(preflight_path.read_text(encoding="utf-8"))
    if not preflight.get("passed") or preflight.get("research_code_revision") != revision:
        raise Study3AnalysisError("preflight did not pass on the E09 execution revision")
    audit = {
        "cells": len(rows),
        "tasks": len(tasks),
        "unique_run_ids": len(run_ids),
        "unique_execution_identities": len(identities),
        "execution_revision": revision,
        "model_keys": dict(sorted(model_keys.items())),
        "model_responses": responses,
        "cells_per_model": dict(sorted(Counter(r["model_id"] for r in rows).items())),
        "cells_per_interface": dict(sorted(Counter(r["interface_id"] for r in rows).items())),
        "raw_manifest_and_metrics_sha256": raw_digest.hexdigest(),
        "raw_manifest_and_metrics_file_count": 1080,
        "infrastructure_attempts": len(
            list((root / "results/infrastructure_attempts/E09").rglob("archive_record.json"))
        ),
    }
    return rows, audit


def summarize(rows: Sequence[dict[str, Any]]) -> list[dict[str, Any]]:
    rng = random.Random(BOOTSTRAP_SEED)
    result: list[dict[str, Any]] = []
    for model in MODELS:
        for interface in INTERFACES:
            group = [r for r in rows if r["model_id"] == model and r["interface_id"] == interface]
            accepted = [float(r["accepted_edit_cell"]) for r in group]
            low, high = bootstrap_mean_ci(accepted, rng)
            attempts = sum(int(r["edit_attempts"]) for r in group)
            acceptances = sum(int(r["edit_acceptances"]) for r in group)
            result.append(
                {
                    "model_id": model,
                    "interface_id": interface,
                    "n": len(group),
                    "accepted_edit_count": int(sum(accepted)),
                    "accepted_edit_rate": statistics.fmean(accepted),
                    "accepted_ci_low": low,
                    "accepted_ci_high": high,
                    "edit_attempts": attempts,
                    "edit_acceptances": acceptances,
                    "attempt_acceptance_rate": acceptances / attempts if attempts else 0.0,
                    "applicable_count": sum(int(r["applicable_final_patch"]) for r in group),
                    "applicable_rate": statistics.fmean(float(r["applicable_final_patch"]) for r in group),
                    "resolved_count": sum(int(r["resolved_at_1"]) for r in group),
                    "resolved_rate": statistics.fmean(float(r["resolved_at_1"]) for r in group),
                    "exact_file_rate": statistics.fmean(float(r["exact_modified_file_match"]) for r in group),
                    "mean_total_tokens": statistics.fmean(float(r["total_tokens"]) for r in group),
                    "mean_elapsed_seconds": statistics.fmean(float(r["elapsed_seconds"]) for r in group),
                    "mean_model_calls": statistics.fmean(float(r["model_calls"]) for r in group),
                    "mean_tool_calls": statistics.fmean(float(r["tool_calls"]) for r in group),
                    "protocol_violation_rate": statistics.fmean(float(r["protocol_violation_count"] > 0) for r in group),
                }
            )
    return result


def _paired(
    rows: Sequence[dict[str, Any]], name: str, endpoint: str, model: str, left: str, right: str,
    family: str,
) -> dict[str, Any]:
    by = {(r["task_id"], r["model_id"], r["interface_id"]): r for r in rows}
    tasks = sorted({r["task_id"] for r in rows})
    left_values = [int(by[(task, model, left)][endpoint]) for task in tasks]
    right_values = [int(by[(task, model, right)][endpoint]) for task in tasks]
    differences = [a - b for a, b in zip(left_values, right_values)]
    rng = random.Random(BOOTSTRAP_SEED + sum(ord(c) for c in name))
    low, high = bootstrap_mean_ci(differences, rng)
    n10, n01, p_value = exact_mcnemar(left_values, right_values)
    return {
        "contrast": name,
        "family": family,
        "endpoint": endpoint,
        "model_id": model,
        "left_interface": left,
        "right_interface": right,
        "tasks": len(tasks),
        "left_count": sum(left_values),
        "right_count": sum(right_values),
        "paired_risk_difference": statistics.fmean(differences),
        "risk_difference_ci_low": low,
        "risk_difference_ci_high": high,
        "discordant_left_only": n10,
        "discordant_right_only": n01,
        "mcnemar_p": p_value,
        "mcnemar_p_holm": p_value,
    }


def contrasts(rows: Sequence[dict[str, Any]]) -> list[dict[str, Any]]:
    primary = _paired(rows, *PRIMARY, family="confirmatory_primary")
    secondary: list[dict[str, Any]] = []
    pairs = (("P002", "P001"), ("P003", "P001"), ("P003", "P002"))
    index = 1
    for endpoint in ("accepted_edit_cell", "resolved_at_1"):
        for model in MODELS:
            for left, right in pairs:
                if endpoint == "accepted_edit_cell" and model == "M003" and (left, right) == (
                    "P002",
                    "P001",
                ):
                    continue
                secondary.append(
                    _paired(
                        rows,
                        f"S{index:02d}_{model}_{left}_vs_{right}_{endpoint}",
                        endpoint,
                        model,
                        left,
                        right,
                        family="prespecified_secondary",
                    )
                )
                index += 1
    adjusted = holm_adjust([row["mcnemar_p"] for row in secondary])
    for row, value in zip(secondary, adjusted):
        row["mcnemar_p_holm"] = value
    return [primary, *secondary]


def compatibility_gate(summaries: Sequence[dict[str, Any]]) -> dict[str, Any]:
    result: dict[str, Any] = {
        "criteria": {
            "minimum_accepted_edit_cells": 12,
            "minimum_attempt_acceptance_rate": 0.5,
            "selection_order": [
                "accepted_edit_rate descending",
                "applicable_rate descending",
                "mean_total_tokens ascending",
                "interface_id ascending",
            ],
            "hidden_test_resolution_used": False,
        },
        "models": {},
    }
    for model in MODELS:
        values = [row for row in summaries if row["model_id"] == model]
        qualifying = [
            row
            for row in values
            if row["accepted_edit_count"] >= 12 and row["attempt_acceptance_rate"] >= 0.5
        ]
        qualifying.sort(
            key=lambda row: (
                -row["accepted_edit_rate"],
                -row["applicable_rate"],
                row["mean_total_tokens"],
                row["interface_id"],
            )
        )
        result["models"][model] = {
            "qualifying_interfaces": [row["interface_id"] for row in qualifying],
            "selected_interface": qualifying[0]["interface_id"] if qualifying else None,
            "eligible_for_fresh_retrieval": bool(qualifying),
        }
    return result


def hierarchical_model(rows: Sequence[dict[str, Any]]) -> dict[str, Any]:
    try:
        from statsmodels.genmod.bayes_mixed_glm import BinomialBayesMixedGLM

        frame = pd.DataFrame(rows)
        formula = (
            "accepted_edit_cell ~ C(interface_id, Treatment(reference='P001')) * "
            "C(model_id, Treatment(reference='M002')) + "
            "C(repository_id, Treatment(reference='R001'))"
        )
        fitted = BinomialBayesMixedGLM.from_formula(
            formula, {"task": "0 + C(task_id)"}, frame
        ).fit_vb()
        names = list(fitted.model.exog_names)
        means = np.asarray(fitted.params[: len(names)], dtype=float)
        standard = np.asarray(fitted.fe_sd, dtype=float)
        coefficients = [
            {
                "term": name,
                "log_odds_mean": float(mean),
                "log_odds_sd": float(sd),
                "odds_ratio": float(math.exp(mean)),
                "or_interval_low": float(math.exp(mean - 1.96 * sd)),
                "or_interval_high": float(math.exp(mean + 1.96 * sd)),
            }
            for name, mean, sd in zip(names, means, standard)
        ]
        return {
            "status": "converged",
            "model_type": "Bayesian binomial mixed model (variational Bayes)",
            "formula": formula,
            "task_random_intercept": True,
            "coefficients": coefficients,
        }
    except Exception as exc:
        return {
            "status": "failed",
            "model_type": "Bayesian binomial mixed model (variational Bayes)",
            "error": repr(exc),
            "coefficients": [],
        }


def failure_analysis(rows: Sequence[dict[str, Any]]) -> dict[str, Any]:
    return {
        "failure_stages": dict(sorted(Counter(r["failure_stage"] for r in rows).items())),
        "finished_reasons": dict(sorted(Counter(r["finished_reason"] for r in rows).items())),
        "cells_with_edit_attempt": sum(r["edit_attempts"] > 0 for r in rows),
        "accepted_edit_cells": sum(r["accepted_edit_cell"] for r in rows),
        "applicable_final_patches": sum(r["applicable_final_patch"] for r in rows),
        "resolved_cells": sum(r["resolved_at_1"] for r in rows),
        "protocol_violation_cells": sum(r["protocol_violation_count"] > 0 for r in rows),
    }


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


def paper_tables(
    path: Path,
    summaries: Sequence[dict[str, Any]],
    contrast_rows: Sequence[dict[str, Any]],
    gate: dict[str, Any],
) -> None:
    lines = [
        "% Generated by scripts/analyze_study3.py; do not edit.",
        r"\begin{table*}[t]",
        r"\centering\small",
        r"\caption{E09 protocol outcomes by model and edit interface.}",
        r"\label{tab:e09-main}",
        r"\begin{tabular}{lllrrrrrr}",
        r"\toprule",
        r"Model & Interface & Accepted & Attempt acceptance & Applicable & Resolved & Exact files & Mean tokens \\",
        r"\midrule",
    ]
    for row in summaries:
        lines.append(
            f"{row['model_id']} & {row['interface_id']} & "
            f"{row['accepted_edit_count']}/{row['n']} ({row['accepted_edit_rate']:.3f}) & "
            f"{row['edit_acceptances']}/{row['edit_attempts']} ({row['attempt_acceptance_rate']:.3f}) & "
            f"{row['applicable_count']}/{row['n']} & {row['resolved_count']}/{row['n']} & "
            f"{row['exact_file_rate']:.3f} & {row['mean_total_tokens']:.0f} \\\\"
        )
    lines.extend([r"\bottomrule", r"\end{tabular}", r"\end{table*}", ""])
    primary = contrast_rows[0]
    lines.extend(
        [
            r"\begin{table}[t]",
            r"\centering\small",
            r"\caption{E09 sole confirmatory paired contrast.}",
            r"\label{tab:e09-primary}",
            r"\begin{tabular}{lr}",
            r"\toprule",
            r"Statistic & Value \\",
            r"\midrule",
            f"P002/P001 accepted & {primary['left_count']}/60 vs {primary['right_count']}/60 \\\\ ",
            f"Paired risk difference & {primary['paired_risk_difference']:+.3f} "
            f"[{primary['risk_difference_ci_low']:+.3f}, {primary['risk_difference_ci_high']:+.3f}] \\\\ ",
            f"Discordant P002/P001 only & {primary['discordant_left_only']}/{primary['discordant_right_only']} \\\\ ",
            f"Exact McNemar $p$ & {primary['mcnemar_p']:.4g} \\\\ ",
            r"\bottomrule",
            r"\end{tabular}",
            r"\end{table}",
            "",
            "% Compatibility gate: " + json.dumps(gate, sort_keys=True),
        ]
    )
    path.write_text("\n".join(lines), encoding="utf-8")


def plots(
    output: Path,
    summaries: Sequence[dict[str, Any]],
    contrast_rows: Sequence[dict[str, Any]],
) -> list[Path]:
    x = np.arange(len(INTERFACES))
    width = 0.24
    colors = {"M002": "#2364aa", "M003": "#f18f01", "M004": "#6a994e"}
    fig, axes = plt.subplots(1, 3, figsize=(11.0, 3.8), sharey=True)
    for axis, endpoint, title in zip(
        axes,
        ("accepted_edit_rate", "applicable_rate", "resolved_rate"),
        ("Executor-accepted edit", "Applicable final patch", "Resolved@1"),
    ):
        for offset_index, model in enumerate(MODELS):
            values = [
                next(
                    row[endpoint]
                    for row in summaries
                    if row["model_id"] == model and row["interface_id"] == interface
                )
                for interface in INTERFACES
            ]
            offset = (offset_index - 1) * width
            axis.bar(x + offset, values, width, label=model, color=colors[model])
        axis.set_xticks(x, INTERFACES)
        axis.set_title(title)
        axis.set_ylim(0, 1)
        axis.grid(axis="y", alpha=0.25)
    axes[0].set_ylabel("Cell rate")
    axes[-1].legend(frameon=False)
    fig.tight_layout()
    bars = output / "e09_protocol_outcomes.pdf"
    fig.savefig(bars)
    fig.savefig(output / "e09_protocol_outcomes.png", dpi=240)
    plt.close(fig)

    accepted = [row for row in contrast_rows if row["endpoint"] == "accepted_edit_cell"]
    forest = list(reversed(accepted))
    estimates = np.asarray([row["paired_risk_difference"] for row in forest])
    low = np.asarray([row["risk_difference_ci_low"] for row in forest])
    high = np.asarray([row["risk_difference_ci_high"] for row in forest])
    fig, axis = plt.subplots(figsize=(8.0, 5.2))
    y = np.arange(len(forest))
    axis.errorbar(
        estimates,
        y,
        xerr=[estimates - low, high - estimates],
        fmt="o",
        color="#2a6f97",
        capsize=3,
    )
    axis.axvline(0, color="black", linewidth=0.8)
    axis.set_yticks(y, [row["contrast"] for row in forest], fontsize=8)
    axis.set_xlabel("Paired accepted-edit risk difference")
    axis.set_xlim(-1, 1)
    axis.grid(axis="x", alpha=0.25)
    fig.tight_layout()
    forest_path = output / "e09_accepted_contrast_forest.pdf"
    fig.savefig(forest_path)
    fig.savefig(output / "e09_accepted_contrast_forest.png", dpi=240)
    plt.close(fig)
    return [bars, forest_path]


def analyze(root: Path) -> dict[str, Any]:
    rows, audit = discover(root)
    summaries = summarize(rows)
    contrast_rows = contrasts(rows)
    gate = compatibility_gate(summaries)
    hierarchy = hierarchical_model(rows)
    failures = failure_analysis(rows)
    output = root / "results/derived/study3"
    output.mkdir(parents=True, exist_ok=True)
    csv_files: list[Path] = []
    for name, values in (
        ("e09_cells.csv", rows),
        ("e09_model_interface_summary.csv", summaries),
        ("e09_contrasts.csv", contrast_rows),
        ("e09_hierarchical_coefficients.csv", hierarchy["coefficients"]),
    ):
        path = output / name
        write_csv(path, values)
        csv_files.append(path)
    table_path = output / "e09_paper_tables.tex"
    paper_tables(table_path, summaries, contrast_rows, gate)
    figures = plots(output, summaries, contrast_rows)
    analysis_revision = git_output(root, "rev-parse", "HEAD")
    manifest = {
        "schema_version": 1,
        "experiment_id": "E09",
        "execution_revision": audit["execution_revision"],
        "analysis_code_revision": analysis_revision,
        "analysis_script_sha256": sha256_file(root / "scripts/analyze_study3.py"),
        "input_cells": 540,
        "input_run_ids": sorted(row["run_id"] for row in rows),
        "raw_manifest_and_metrics_sha256": audit["raw_manifest_and_metrics_sha256"],
        "bootstrap_samples": BOOTSTRAPS,
        "bootstrap_seed": BOOTSTRAP_SEED,
        "primary_test": "P002 versus P001 accepted-edit cell on M003; two-sided exact McNemar",
        "secondary_multiplicity": "Holm across 17 prespecified binary contrasts",
        "gate_uses_resolution": False,
    }
    report = {
        "schema_version": 1,
        "experiment_id": "E09",
        "claim_boundary": {
            "confirmatory": "P002 versus P001 accepted-edit cell on M003 only",
            "resolution": "secondary end-to-end endpoint",
            "task_reuse": "same 60 E08 tasks; new protocol-causality estimand, not independent benchmark replication",
            "language_limit": "repository and language remain confounded",
        },
        "audit": audit,
        "model_interface_summaries": summaries,
        "contrasts": contrast_rows,
        "compatibility_gate": gate,
        "hierarchical_model": hierarchy,
        "failure_analysis": failures,
        "analysis_manifest": manifest,
    }
    manifest_path = output / "analysis_manifest.json"
    report_path = output / "e09_analysis.json"
    manifest_path.write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n")
    report_path.write_text(json.dumps(report, indent=2, sort_keys=True) + "\n")
    generated = [*csv_files, table_path, *figures, manifest_path, report_path]
    checksums = {path.name: sha256_file(path) for path in generated}
    (output / "SHA256SUMS.json").write_text(
        json.dumps(checksums, indent=2, sort_keys=True) + "\n"
    )
    return {
        **report,
        "output_directory": str(output),
        "generated_file_count": len(generated) + 1,
        "checksums": checksums,
    }


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--root", type=Path, default=Path(__file__).resolve().parents[1])
    arguments = parser.parse_args()
    try:
        result = analyze(arguments.root.resolve())
    except (
        Study3AnalysisError,
        OSError,
        ValueError,
        KeyError,
        subprocess.CalledProcessError,
    ) as exc:
        print(f"STUDY 3 ANALYSIS FAILED: {exc}")
        return 1
    print(json.dumps(result, indent=2, sort_keys=True))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())