scripts/analyze_study2.py

#!/usr/bin/env python3
"""Prespecified analysis for the frozen E08 Study 2 execution.

The script refuses incomplete, duplicated, mixed-revision, or off-profile data.
It keeps the temperature-zero main study separate from the temperature-0.2
three-seed reliability extension and emits the paper's machine-generated
tables, figures, checksums, and analysis manifest.
"""

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
import subprocess
import tomllib
from typing import Any, Iterable, Sequence

import matplotlib

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


TREATMENTS = ("H000", "H003", "H007", "H011", "H018", "A001", "A002")
NON_ORACLE_TREATMENTS = ("H000", "H003", "H007", "H011", "A001", "A002")
MODELS = ("M002", "M003")
REPOSITORIES = ("R001", "R002", "R003")
EXPECTED_REVISION = "58933d6fa8af09fcc5a832fb3b523ffb4182bc50"
MAIN_MODEL_HASHES = {
    "M002": "8b11f5093d4eb05ddf7019f7edd2500bcda76ba4c87646690117c117a2d62240",
    "M003": "6001206e5358792f197c1edaa21ebccd2b3163931224aee1beec460378a8fb10",
}
RELIABILITY_MODEL_HASHES = {
    "M002": "ac50d60bd4095b23bfe8f55ff76144c608fd4ebf81fa601c7c97576c9a8f3204",
    "M003": "28596cfc56fd1b18f3fbb0d8f49e03129c03ca3302775d06e2e206d8d273eeb3",
}
BOOTSTRAPS = 20_000
BOOTSTRAP_SEED = 20260718

# The confirmatory contrast is deliberately outside the secondary Holm family.
PRIMARY_CONTRAST = ("P1_H007_vs_H000_M002", "M002", "H007", "H000")
SECONDARY_CONTRASTS = (
    ("S01_H007_vs_H000_M003", "M003", "H007", "H000"),
    ("S02_H003_vs_H000_M002", "M002", "H003", "H000"),
    ("S03_H003_vs_H000_M003", "M003", "H003", "H000"),
    ("S04_H007_vs_H003_M002", "M002", "H007", "H003"),
    ("S05_H007_vs_H003_M003", "M003", "H007", "H003"),
    ("S06_H011_vs_H007_M002", "M002", "H011", "H007"),
    ("S07_H011_vs_H007_M003", "M003", "H011", "H007"),
    ("S08_H018_vs_H007_M002", "M002", "H018", "H007"),
    ("S09_H018_vs_H007_M003", "M003", "H018", "H007"),
    ("S10_A001_vs_H007_M002", "M002", "A001", "H007"),
    ("S11_A001_vs_H007_M003", "M003", "A001", "H007"),
    ("S12_A002_vs_H007_M002", "M002", "A002", "H007"),
    ("S13_A002_vs_H007_M003", "M003", "A002", "H007"),
)


class Study2AnalysisError(RuntimeError):
    """Raised when the frozen analysis contract cannot be satisfied."""


def percentile(values: Sequence[float], probability: float) -> float:
    if not values:
        raise ValueError("percentile requires at least one value")
    ordered = sorted(float(value) for value in 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]:
    """Percentile CI from resampling the task-level values."""
    if not values:
        return math.nan, math.nan
    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]:
    """Two-sided exact McNemar/binomial test for matched binary outcomes."""
    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.0 * tail)


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


def fleiss_kappa(matrix: Sequence[Sequence[int]]) -> float:
    """Fleiss' kappa for binary counts [failures, successes] per cell."""
    values = np.asarray(matrix, dtype=float)
    if values.ndim != 2 or values.shape[1] != 2:
        raise ValueError("Fleiss matrix must be N x 2")
    raters = values.sum(axis=1)
    if not np.all(raters == raters[0]) or raters[0] < 2:
        raise ValueError("Fleiss rows must have a constant rater count >=2")
    n = float(raters[0])
    observed = np.mean((np.square(values).sum(axis=1) - n) / (n * (n - 1.0)))
    marginal = values.sum(axis=0) / values.sum()
    expected = float(np.square(marginal).sum())
    if math.isclose(expected, 1.0):
        return math.nan
    return (float(observed) - expected) / (1.0 - expected)


def krippendorff_alpha_nominal(matrix: Sequence[Sequence[int]]) -> float:
    """Krippendorff's alpha for complete binary nominal ratings."""
    values = np.asarray(matrix, dtype=int)
    if values.ndim != 2:
        raise ValueError("ratings must be a two-dimensional matrix")
    total_pairs = values.shape[0] * values.shape[1] * (values.shape[1] - 1)
    disagree = 0
    for row in values:
        zeros = int(np.sum(row == 0))
        ones = int(np.sum(row == 1))
        disagree += 2 * zeros * ones
    observed = disagree / total_pairs
    flat = values.ravel()
    zeros = int(np.sum(flat == 0))
    ones = int(np.sum(flat == 1))
    expected = 2 * zeros * ones / (len(flat) * (len(flat) - 1))
    return math.nan if expected == 0 else 1.0 - observed / expected


def icc_one_way(matrix: Sequence[Sequence[int]]) -> float:
    """One-way random-effects single-measure ICC(1,1)."""
    values = np.asarray(matrix, dtype=float)
    groups, raters = values.shape
    if groups < 2 or raters < 2:
        raise ValueError("ICC requires at least two groups and two ratings")
    means = values.mean(axis=1)
    ms_between = raters * float(np.var(means, ddof=1))
    ms_within = float(np.square(values - means[:, None]).sum() / (groups * (raters - 1)))
    denominator = ms_between + (raters - 1) * ms_within
    return math.nan if denominator == 0 else (ms_between - ms_within) / denominator


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


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


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


def load_task_metadata(root: Path) -> dict[str, dict[str, Any]]:
    repositories: dict[str, dict[str, Any]] = {}
    url_to_repo: dict[str, str] = {}
    for path in sorted((root / "configs/repositories").glob("R*.toml")):
        data = tomllib.loads(path.read_text(encoding="utf-8"))
        repositories[str(data["repository_id"])] = data
        url_to_repo[str(data["repository_url"])] = str(data["repository_id"])
    tasks: dict[str, dict[str, Any]] = {}
    frozen = {
        line.strip()
        for line in (root / "tasks/splits/study2_confirmatory.txt").read_text(encoding="utf-8").splitlines()
        if line.strip() and not line.lstrip().startswith("#")
    }
    for task_id in sorted(frozen):
        path = root / "tasks/manifests" / f"{task_id}.toml"
        data = tomllib.loads(path.read_text(encoding="utf-8"))
        repository_id = url_to_repo.get(str(data["repository_url"]))
        if repository_id is None:
            raise Study2AnalysisError(f"task {task_id} has an unknown repository URL")
        tasks[task_id] = {
            **data,
            "repository_id": repository_id,
            "repository_name": repositories[repository_id]["name"],
        }
    if len(tasks) != 60 or Counter(x["repository_id"] for x in tasks.values()) != Counter(
        {"R001": 20, "R002": 20, "R003": 20}
    ):
        raise Study2AnalysisError("frozen task metadata is not the registered 20/20/20 split")
    return tasks


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


def discover(root: Path) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    """Load E08 results and enforce the complete frozen execution identity."""
    tasks = load_task_metadata(root)
    paths = sorted((root / "results/raw/E08").rglob("final_metrics.json"))
    rows: list[dict[str, Any]] = []
    identities: set[tuple[Any, ...]] = set()
    run_ids: set[str] = set()
    raw_digest = sha256()
    model_keys: Counter[str] = Counter()
    model_hashes: Counter[str] = Counter()
    revisions: Counter[str] = Counter()
    response_count = 0
    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"]
        task_id = str(identity["task_id"])
        if task_id not in tasks:
            raise Study2AnalysisError(f"unexpected task {task_id}")
        key = tuple(
            identity[name]
            for name in ("task_id", "harness_id", "model_id", "seed", "repetition")
        )
        if key in identities:
            raise Study2AnalysisError(f"duplicate execution identity {key}")
        identities.add(key)
        run_id = str(final["run_id"])
        if run_id in run_ids:
            raise Study2AnalysisError(f"duplicate run ID {run_id}")
        run_ids.add(run_id)
        if final["task_id"] != task_id or final["harness_id"] != identity["harness_id"]:
            raise Study2AnalysisError(f"manifest/final identity mismatch in {directory}")
        model = manifest["resolved_model"]["agent_model"]
        runtime = manifest["resolved_model"]["agent_runtime"]
        native = runtime["native_record"]
        if runtime["inference_key"] != identity["model_key"]:
            raise Study2AnalysisError(f"runtime/model-key mismatch in {directory}")
        loaded = native.get("loaded_instances") or []
        if len(loaded) != 1 or int(loaded[0]["config"]["context_length"]) != 65536:
            raise Study2AnalysisError(f"wrong runtime residency/context in {directory}")
        if any(len(item.get("after_instances", [])) != 1 for item in final["residency_transitions"]):
            raise Study2AnalysisError(f"non-exclusive model residency in {directory}")
        responses = sorted(directory.glob("model_response_*.json"))
        if len(responses) != int(final["model_calls"]):
            raise Study2AnalysisError(f"model response count mismatch in {directory}")
        response_count += len(responses)
        revisions[str(identity["code_revision"])] += 1
        model_keys[str(identity["model_key"])] += 1
        model_hashes[str(identity["model_config_hash"])] += 1
        gold_files = set(str(value) for value in tasks[task_id]["gold_files"])
        modified_files = set(str(value) for value in final["modified_files"])
        usage = final.get("usage") or {}
        row = {
            "run_id": run_id,
            "task_id": task_id,
            "repository_id": tasks[task_id]["repository_id"],
            "repository_name": tasks[task_id]["repository_name"],
            "language": tasks[task_id]["language"],
            "difficulty": tasks[task_id]["difficulty"],
            "treatment_id": str(identity["harness_id"]),
            "model_id": str(identity["model_id"]),
            "model_key": str(identity["model_key"]),
            "seed": int(identity["seed"]),
            "repetition": int(identity["repetition"]),
            "temperature": float(model["temperature"]),
            "top_p": float(model["top_p"]),
            "model_config_hash": str(identity["model_config_hash"]),
            "code_revision": str(identity["code_revision"]),
            "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"])),
            "patch_applied": int(bool(final["patch_applied"])),
            "empty_patch": int(not modified_files),
            "exact_modified_file_match": int(modified_files == gold_files),
            "gold_file_modified_recall": len(modified_files & gold_files) / len(gold_files),
            "search_file_recall_at_10": float(final["search_localization_metrics"]["file_recall_at_10"]),
            "read_file_recall_at_10": float(final["read_localization_metrics"]["file_recall_at_10"]),
            "search_all_gold": int(bool(final["search_localization_metrics"]["all_gold_in_top_10"])),
            "read_all_gold": int(bool(final["read_localization_metrics"]["all_gold_in_top_10"])),
            "model_calls": int(final["model_calls"]),
            "tool_calls": int(final["tool_calls"]),
            "test_runs": int(final["test_runs"]),
            "elapsed_seconds": float(final["elapsed_seconds"]),
            "model_elapsed_seconds": float(final["model_elapsed_seconds"]),
            "prompt_tokens": int(usage["prompt_tokens"]),
            "completion_tokens": int(usage["completion_tokens"]),
            "total_tokens": int(usage["total_tokens"]),
            "model_switch_count": int(final["model_switch_count"]),
            "model_switch_seconds": float(final["model_switch_seconds"]),
            "protocol_violation_count": len(final["protocol_violations"]),
            "protocol_violations": list(final["protocol_violations"]),
            "finished_reason": str(final["finished_reason"]),
            "failure_stage": str(final["failure_stage"]),
            "gold_file_count": len(gold_files),
            "modified_file_count": len(modified_files),
        }
        rows.append(row)
        for artifact in (directory / "run_manifest.json", path):
            relative = artifact.relative_to(root).as_posix().encode()
            raw_digest.update(relative + b"\0" + artifact.read_bytes() + b"\0")

    main = [row for row in rows if row["repetition"] == 0]
    reliability = [row for row in rows if row["repetition"] == 1]
    expected_main = {
        (task, treatment, model, 0, 0)
        for task in tasks
        for treatment in TREATMENTS
        for model in MODELS
    }
    observed_main = {
        (r["task_id"], r["treatment_id"], r["model_id"], r["seed"], r["repetition"])
        for r in main
    }
    if len(rows) != 912 or len(main) != 840 or len(reliability) != 72:
        raise Study2AnalysisError(
            f"E08 count mismatch: all={len(rows)}, main={len(main)}, reliability={len(reliability)}"
        )
    if observed_main != expected_main:
        raise Study2AnalysisError(f"main grid mismatch; missing={sorted(expected_main-observed_main)[:5]}")
    repeat_manifest = json.loads(
        (root / "configs/reliability/E08_repeat_cells.json").read_text(encoding="utf-8")
    )
    expected_reliability = {
        (cell["task_id"], cell["treatment_id"], cell["model_id"], seed, 1)
        for cell in repeat_manifest["cells"]
        for seed in repeat_manifest["seeds"]
    }
    observed_reliability = {
        (r["task_id"], r["treatment_id"], r["model_id"], r["seed"], r["repetition"])
        for r in reliability
    }
    if observed_reliability != expected_reliability:
        raise Study2AnalysisError("reliability grid differs from its frozen manifest")
    if revisions != Counter({EXPECTED_REVISION: 912}):
        raise Study2AnalysisError(f"mixed/wrong execution revisions: {dict(revisions)}")
    if any(r["temperature"] != 0.0 or r["top_p"] != 1.0 for r in main):
        raise Study2AnalysisError("main generation profile differs from preregistration")
    if any(r["temperature"] != 0.2 or r["top_p"] != 1.0 for r in reliability):
        raise Study2AnalysisError("reliability generation profile differs from amendment PA-005")
    for model_id, expected_hash in MAIN_MODEL_HASHES.items():
        if {r["model_config_hash"] for r in main if r["model_id"] == model_id} != {expected_hash}:
            raise Study2AnalysisError(f"wrong main model hash for {model_id}")
    for model_id, expected_hash in RELIABILITY_MODEL_HASHES.items():
        if {r["model_config_hash"] for r in reliability if r["model_id"] == model_id} != {expected_hash}:
            raise Study2AnalysisError(f"wrong reliability model hash for {model_id}")

    audit = {
        "cells": len(rows),
        "main_cells": len(main),
        "reliability_cells": len(reliability),
        "unique_run_ids": len(run_ids),
        "unique_execution_identities": len(identities),
        "execution_revision": EXPECTED_REVISION,
        "model_keys": dict(sorted(model_keys.items())),
        "model_config_hashes": dict(sorted(model_hashes.items())),
        "model_responses": response_count,
        "task_count": len(tasks),
        "main_cells_per_task": dict(Counter(Counter(r["task_id"] for r in main).values())),
        "reliability_groups": len(
            {(r["task_id"], r["treatment_id"], r["model_id"]) for r in reliability}
        ),
        "reliability_cells_per_group": dict(
            Counter(
                Counter(
                    (r["task_id"], r["treatment_id"], r["model_id"])
                    for r in reliability
                ).values()
            )
        ),
        "infrastructure_failure_archives": len(
            list((root / "results/raw/E08").glob("**/attempts/infrastructure/*.json"))
        ),
        "raw_manifest_and_metrics_sha256": raw_digest.hexdigest(),
        "raw_manifest_and_metrics_file_count": len(rows) * 2,
    }
    if audit["infrastructure_failure_archives"] != 0:
        raise Study2AnalysisError("unexpected infrastructure failure archive")
    return rows, audit


def _mean(group: Sequence[dict[str, Any]], key: str) -> float:
    return statistics.fmean(float(row[key]) for row in group)


def _median(group: Sequence[dict[str, Any]], key: str) -> float:
    return statistics.median(float(row[key]) for row in group)


def summarize_main(main: Sequence[dict[str, Any]]) -> list[dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    rng = random.Random(BOOTSTRAP_SEED)
    for treatment in TREATMENTS:
        for model in MODELS:
            group = [r for r in main if r["treatment_id"] == treatment and r["model_id"] == model]
            resolved = [float(r["resolved_at_1"]) for r in group]
            low, high = bootstrap_mean_ci(resolved, rng)
            rows.append(
                {
                    "treatment_id": treatment,
                    "model_id": model,
                    "n": len(group),
                    "resolved_count": int(sum(resolved)),
                    "resolved_rate": statistics.fmean(resolved),
                    "resolved_ci_low": low,
                    "resolved_ci_high": high,
                    "fail_to_pass_rate": _mean(group, "fail_to_pass"),
                    "pass_to_pass_rate": _mean(group, "pass_to_pass"),
                    "patch_apply_rate": _mean(group, "patch_applied"),
                    "empty_patch_rate": _mean(group, "empty_patch"),
                    "exact_modified_file_match_rate": _mean(group, "exact_modified_file_match"),
                    "gold_file_modified_recall": _mean(group, "gold_file_modified_recall"),
                    "search_file_recall_at_10": _mean(group, "search_file_recall_at_10"),
                    "read_file_recall_at_10": _mean(group, "read_file_recall_at_10"),
                    "mean_model_calls": _mean(group, "model_calls"),
                    "mean_tool_calls": _mean(group, "tool_calls"),
                    "mean_test_runs": _mean(group, "test_runs"),
                    "mean_prompt_tokens": _mean(group, "prompt_tokens"),
                    "mean_completion_tokens": _mean(group, "completion_tokens"),
                    "mean_total_tokens": _mean(group, "total_tokens"),
                    "median_elapsed_seconds": _median(group, "elapsed_seconds"),
                    "mean_elapsed_seconds": _mean(group, "elapsed_seconds"),
                    "mean_model_switch_count": _mean(group, "model_switch_count"),
                    "mean_model_switch_seconds": _mean(group, "model_switch_seconds"),
                    "protocol_violation_rate": statistics.fmean(
                        float(r["protocol_violation_count"] > 0) for r in group
                    ),
                }
            )
    return rows


def repository_strata(main: Sequence[dict[str, Any]]) -> list[dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    for repository in REPOSITORIES:
        for treatment in TREATMENTS:
            for model in MODELS:
                group = [
                    r
                    for r in main
                    if r["repository_id"] == repository
                    and r["treatment_id"] == treatment
                    and r["model_id"] == model
                ]
                rows.append(
                    {
                        "repository_id": repository,
                        "language": group[0]["language"],
                        "treatment_id": treatment,
                        "model_id": model,
                        "n": len(group),
                        "resolved_count": sum(r["resolved_at_1"] for r in group),
                        "resolved_rate": _mean(group, "resolved_at_1"),
                        "patch_apply_rate": _mean(group, "patch_applied"),
                        "search_file_recall_at_10": _mean(group, "search_file_recall_at_10"),
                        "read_file_recall_at_10": _mean(group, "read_file_recall_at_10"),
                        "mean_total_tokens": _mean(group, "total_tokens"),
                        "mean_elapsed_seconds": _mean(group, "elapsed_seconds"),
                    }
                )
    return rows


def paired_contrast(
    main: Sequence[dict[str, Any]],
    name: str,
    model: str,
    left_treatment: str,
    right_treatment: str,
    family: str,
    rng: random.Random,
) -> dict[str, Any]:
    lookup = {
        (r["task_id"], r["treatment_id"]): r for r in main if r["model_id"] == model
    }
    tasks = sorted({r["task_id"] for r in main if r["model_id"] == model})
    left = [lookup[(task, left_treatment)] for task in tasks]
    right = [lookup[(task, right_treatment)] for task in tasks]
    a = [r["resolved_at_1"] for r in left]
    b = [r["resolved_at_1"] for r in right]
    differences = [float(x - y) for x, y in zip(a, b)]
    low, high = bootstrap_mean_ci(differences, rng)
    n10, n01, p_value = exact_mcnemar(a, b)
    result: dict[str, Any] = {
        "contrast": name,
        "family": family,
        "model_id": model,
        "left_treatment": left_treatment,
        "right_treatment": right_treatment,
        "tasks": len(tasks),
        "left_resolved": sum(a),
        "right_resolved": sum(b),
        "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": math.nan if family == "secondary" else p_value,
    }
    for metric in ("total_tokens", "elapsed_seconds", "model_calls", "tool_calls"):
        delta = [float(x[metric] - y[metric]) for x, y in zip(left, right)]
        delta_low, delta_high = bootstrap_mean_ci(delta, rng)
        result[f"mean_{metric}_difference"] = statistics.fmean(delta)
        result[f"{metric}_difference_ci_low"] = delta_low
        result[f"{metric}_difference_ci_high"] = delta_high
    return result


def contrasts(main: Sequence[dict[str, Any]]) -> list[dict[str, Any]]:
    rng = random.Random(BOOTSTRAP_SEED + 1)
    primary = paired_contrast(main, *PRIMARY_CONTRAST, family="confirmatory_primary", rng=rng)
    secondary = [
        paired_contrast(main, *spec, family="secondary", rng=rng) for spec in SECONDARY_CONTRASTS
    ]
    adjusted = holm_adjust([float(row["mcnemar_p"]) for row in secondary])
    for row, value in zip(secondary, adjusted):
        row["mcnemar_p_holm"] = value
    return [primary, *secondary]


def hierarchical_model(main: Sequence[dict[str, Any]]) -> dict[str, Any]:
    """Fit the registered task-random-intercept model, with GEE fallback."""
    data = pd.DataFrame(main).rename(
        columns={"treatment_id": "treatment", "model_id": "model", "repository_id": "repository"}
    )
    formula = (
        "resolved_at_1 ~ C(treatment, Treatment(reference='H000')) * "
        "C(model, Treatment(reference='M002')) + "
        "C(repository, Treatment(reference='R001'))"
    )
    note = (
        "Repository fixed effects encode the registered repository/language strata. "
        "Language is not entered separately because Python occurs in only R003 and is collinear."
    )
    try:
        from statsmodels.genmod.bayes_mixed_glm import BinomialBayesMixedGLM

        model = BinomialBayesMixedGLM.from_formula(
            formula, {"task_intercept": "0 + C(task_id)"}, data
        )
        fit = model.fit_vb()
        converged = bool(fit.optim_retvals.get("success", False))
        if not converged:
            raise RuntimeError(str(fit.optim_retvals.get("message", "VB did not converge")))
        coefficients = []
        for name, estimate, sd in zip(model.exog_names, fit.fe_mean, fit.fe_sd):
            coefficients.append(
                {
                    "term": name,
                    "log_odds": float(estimate),
                    "standard_error_or_posterior_sd": float(sd),
                    "odds_ratio": math.exp(float(estimate)),
                    "interval_low": math.exp(float(estimate) - 1.96 * float(sd)),
                    "interval_high": math.exp(float(estimate) + 1.96 * float(sd)),
                    "p_value": math.nan,
                }
            )
        return {
            "model_type": "Bayesian binomial mixed model (variational Bayes)",
            "formula": formula,
            "random_effects": "task random intercept",
            "converged": True,
            "optimizer_message": str(fit.optim_retvals.get("message", "")),
            "coefficient_interval": "normal approximation to 95% posterior interval",
            "note": note,
            "coefficients": coefficients,
        }
    except Exception as mixed_error:
        import statsmodels.api as sm

        try:
            gee = sm.GEE.from_formula(
                formula,
                groups="task_id",
                data=data,
                family=sm.families.Binomial(),
                cov_struct=sm.cov_struct.Exchangeable(),
            ).fit()
            coefficients = []
            for name in gee.params.index:
                estimate = float(gee.params[name])
                error = float(gee.bse[name])
                coefficients.append(
                    {
                        "term": str(name),
                        "log_odds": estimate,
                        "standard_error_or_posterior_sd": error,
                        "odds_ratio": math.exp(estimate),
                        "interval_low": math.exp(estimate - 1.96 * error),
                        "interval_high": math.exp(estimate + 1.96 * error),
                        "p_value": float(gee.pvalues[name]),
                    }
                )
            return {
                "model_type": "task-clustered binomial GEE (preregistered fallback)",
                "formula": formula,
                "random_effects": None,
                "converged": bool(gee.converged),
                "mixed_model_failure": repr(mixed_error),
                "coefficient_interval": "95% robust Wald confidence interval",
                "note": note,
                "coefficients": coefficients,
            }
        except Exception as gee_error:
            return {
                "model_type": "stratified bootstrap only (preregistered terminal fallback)",
                "formula": formula,
                "random_effects": None,
                "converged": False,
                "mixed_model_failure": repr(mixed_error),
                "gee_failure": repr(gee_error),
                "note": note,
                "coefficients": [],
            }


def reliability_analysis(
    reliability: Sequence[dict[str, Any]], main: Sequence[dict[str, Any]]
) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    grouped: dict[tuple[str, str, str], list[dict[str, Any]]] = defaultdict(list)
    for row in reliability:
        grouped[(row["task_id"], row["treatment_id"], row["model_id"])].append(row)
    main_lookup = {
        (r["task_id"], r["treatment_id"], r["model_id"]): r for r in main
    }
    group_rows: list[dict[str, Any]] = []
    rating_matrix: list[list[int]] = []
    for key in sorted(grouped):
        values = sorted(grouped[key], key=lambda row: row["seed"])
        if [r["seed"] for r in values] != [0, 1, 2]:
            raise Study2AnalysisError(f"reliability seeds are not 0/1/2 for {key}")
        ratings = [r["resolved_at_1"] for r in values]
        rating_matrix.append(ratings)
        pair_agreement = statistics.fmean(
            float(ratings[i] == ratings[j]) for i in range(3) for j in range(i + 1, 3)
        )
        main_value = main_lookup[key]["resolved_at_1"]
        majority = int(sum(ratings) >= 2)
        group_rows.append(
            {
                "task_id": key[0],
                "treatment_id": key[1],
                "model_id": key[2],
                "repository_id": values[0]["repository_id"],
                "seed_0_resolved": ratings[0],
                "seed_1_resolved": ratings[1],
                "seed_2_resolved": ratings[2],
                "successes_across_seeds": sum(ratings),
                "success_rate_across_seeds": statistics.fmean(ratings),
                "sample_variance_across_seeds": statistics.variance(ratings),
                "unanimous": int(len(set(ratings)) == 1),
                "pairwise_agreement": pair_agreement,
                "majority_resolved": majority,
                "main_temperature_zero_resolved": main_value,
                "majority_agrees_with_main": int(majority == main_value),
            }
        )
    counts = [[3 - sum(row), sum(row)] for row in rating_matrix]
    flat = [value for row in rating_matrix for value in row]
    by_model: dict[str, Any] = {}
    by_treatment: dict[str, Any] = {}
    for dimension, target in (("model_id", by_model), ("treatment_id", by_treatment)):
        levels = MODELS if dimension == "model_id" else NON_ORACLE_TREATMENTS
        for level in levels:
            subset = [row for row in group_rows if row[dimension] == level]
            target[level] = {
                "groups": len(subset),
                "unanimous_rate": statistics.fmean(float(r["unanimous"]) for r in subset),
                "pairwise_agreement": statistics.fmean(r["pairwise_agreement"] for r in subset),
                "mean_success_rate": statistics.fmean(r["success_rate_across_seeds"] for r in subset),
            }
    kappa = fleiss_kappa(counts)
    alpha = krippendorff_alpha_nominal(rating_matrix)
    icc = icc_one_way(rating_matrix)
    summary = {
        "temperature": 0.2,
        "seeds": [0, 1, 2],
        "groups": len(group_rows),
        "observations": len(flat),
        "successes": sum(flat),
        "success_rate": statistics.fmean(flat),
        "unanimous_groups": sum(r["unanimous"] for r in group_rows),
        "unanimous_rate": statistics.fmean(float(r["unanimous"]) for r in group_rows),
        "mean_pairwise_agreement": statistics.fmean(r["pairwise_agreement"] for r in group_rows),
        "mean_within_group_binary_variance": statistics.fmean(
            r["sample_variance_across_seeds"] for r in group_rows
        ),
        "fleiss_kappa": None if math.isnan(kappa) else kappa,
        "krippendorff_alpha_nominal": None if math.isnan(alpha) else alpha,
        "icc_1_1": None if math.isnan(icc) else icc,
        "chance_corrected_agreement_note": (
            "Fleiss kappa, Krippendorff alpha, and ICC are undefined when all ratings "
            "occupy one outcome category; raw agreement remains descriptive."
        ),
        "majority_agreement_with_temperature_zero": statistics.fmean(
            float(r["majority_agrees_with_main"]) for r in group_rows
        ),
        "temperature_zero_success_rate_on_matched_cells": statistics.fmean(
            r["main_temperature_zero_resolved"] for r in group_rows
        ),
        "temperature_point_two_mean_success_rate_on_matched_cells": statistics.fmean(
            r["success_rate_across_seeds"] for r in group_rows
        ),
        "by_model": by_model,
        "by_treatment": by_treatment,
        "interpretation_guardrail": (
            "The temperature-0.2 repetitions are a separate reliability sensitivity analysis; "
            "they are not pooled with or substituted for the temperature-zero primary endpoint."
        ),
    }
    return group_rows, summary


def failure_taxonomy(main: Sequence[dict[str, Any]]) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    for treatment in TREATMENTS:
        for model in MODELS:
            group = [r for r in main if r["treatment_id"] == treatment and r["model_id"] == model]
            stages = Counter(r["failure_stage"] for r in group)
            reasons = Counter(r["finished_reason"] for r in group)
            for kind, counts in (("failure_stage", stages), ("finished_reason", reasons)):
                for value, count in sorted(counts.items()):
                    rows.append(
                        {
                            "treatment_id": treatment,
                            "model_id": model,
                            "taxonomy": kind,
                            "value": value,
                            "count": count,
                            "rate": count / len(group),
                        }
                    )
    violations: Counter[str] = Counter()
    for row in main:
        violations.update(str(value) for value in row["protocol_violations"])
    summary = {
        "failure_stages": dict(sorted(Counter(r["failure_stage"] for r in main).items())),
        "finished_reasons": dict(sorted(Counter(r["finished_reason"] for r in main).items())),
        "protocol_violation_cells": sum(r["protocol_violation_count"] > 0 for r in main),
        "protocol_violation_events": sum(r["protocol_violation_count"] for r in main),
        "protocol_violation_types": dict(sorted(violations.items())),
        "empty_patch_cells": sum(r["empty_patch"] for r in main),
        "nonempty_patch_cells": sum(not r["empty_patch"] for r in main),
        "applicable_patch_cells": sum(r["patch_applied"] for r in main),
    }
    return rows, summary


def exploratory_tool_protocol_diagnostic(
    root: Path,
) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    """Derive post-outcome tool-protocol diagnostics from trajectory truth.

    This intentionally does not consume ``final_metrics.tool_counts``. A
    post-outcome audit found that field increments failed calls both on entry
    and on the error path, whereas the aggregate ``tool_calls`` field and one
    event per call in ``trajectory.jsonl`` agree exactly.
    """
    buckets: dict[tuple[str, str], Counter[str]] = defaultdict(Counter)
    tool_names: dict[tuple[str, str], Counter[str]] = defaultdict(Counter)
    for path in sorted((root / "results/raw/E08").rglob("run_manifest.json")):
        manifest = json.loads(path.read_text(encoding="utf-8"))
        identity = manifest["identity"]
        if int(identity["repetition"]) != 0:
            continue
        key = (str(identity["model_id"]), str(identity["harness_id"]))
        for line in (path.parent / "trajectory.jsonl").read_text(encoding="utf-8").splitlines():
            event = json.loads(line)
            if event["event_type"] != "tool_call":
                continue
            payload = event["payload"]
            name = str(payload.get("name"))
            is_error = bool(payload.get("is_error"))
            result = payload.get("result")
            buckets[key]["tool_events"] += 1
            buckets[key]["tool_error_events" if is_error else "tool_executed_events"] += 1
            tool_names[key][name] += 1
            if "unknown" in json.dumps(result, sort_keys=True).lower():
                buckets[key]["unknown_tool_events"] += 1
            if name == "apply_patch":
                buckets[key]["patch_attempts"] += 1
                if is_error:
                    buckets[key]["patch_protocol_errors"] += 1
                else:
                    buckets[key]["patch_executor_calls"] += 1
                    if isinstance(result, dict) and bool(result.get("accepted")):
                        buckets[key]["patch_executor_acceptances"] += 1
        for response_path in path.parent.glob("model_response_*.json"):
            response = json.loads(response_path.read_text(encoding="utf-8"))
            message = (response.get("choices") or [{}])[0].get("message") or {}
            buckets[key]["model_responses"] += 1
            if message.get("tool_calls"):
                buckets[key]["responses_with_tool_calls"] += 1
            else:
                buckets[key]["responses_without_tool_calls"] += 1
                if not message.get("content"):
                    buckets[key]["empty_responses_without_tool_calls"] += 1
    rows: list[dict[str, Any]] = []
    count_fields = (
        "model_responses",
        "responses_with_tool_calls",
        "responses_without_tool_calls",
        "empty_responses_without_tool_calls",
        "tool_events",
        "tool_executed_events",
        "tool_error_events",
        "unknown_tool_events",
        "patch_attempts",
        "patch_protocol_errors",
        "patch_executor_calls",
        "patch_executor_acceptances",
    )
    for model in MODELS:
        for treatment in TREATMENTS:
            key = (model, treatment)
            counts = buckets[key]
            rows.append(
                {
                    "model_id": model,
                    "treatment_id": treatment,
                    **{name: int(counts[name]) for name in count_fields},
                    "tool_error_rate": (
                        counts["tool_error_events"] / counts["tool_events"]
                        if counts["tool_events"]
                        else math.nan
                    ),
                    "patch_protocol_error_rate": (
                        counts["patch_protocol_errors"] / counts["patch_attempts"]
                        if counts["patch_attempts"]
                        else math.nan
                    ),
                    "trajectory_tool_names": json.dumps(
                        dict(sorted(tool_names[key].items())), sort_keys=True
                    ),
                }
            )
    by_model: dict[str, Any] = {}
    for model in MODELS:
        group = [row for row in rows if row["model_id"] == model]
        totals = {key: sum(int(row[key]) for row in group) for key in count_fields}
        totals["tool_error_rate"] = totals["tool_error_events"] / totals["tool_events"]
        totals["patch_protocol_error_rate"] = (
            totals["patch_protocol_errors"] / totals["patch_attempts"]
        )
        by_model[model] = totals
    summary = {
        "status": "post_outcome_exploratory_failure_mechanism_analysis",
        "source_of_truth": "one tool_call event per call in trajectory.jsonl",
        "endpoint_impact": "none; primary and secondary endpoint fields are unaffected",
        "telemetry_caveat": (
            "final_metrics.tool_counts double-counts failed tool invocations; this diagnostic "
            "uses trajectory events, while final_metrics.tool_calls remains correct"
        ),
        "by_model": by_model,
    }
    return rows, summary


def repository_characteristics(root: Path) -> list[dict[str, Any]]:
    audit = json.loads((root / "docs/STUDY2_DESIGN_AUDIT.json").read_text(encoding="utf-8"))
    rows = []
    for item in audit["repositories"]:
        rows.append(
            {
                "repository_id": item["repository_id"],
                "name": item["name"],
                "language": item["language"],
                "tasks": audit["repository_task_counts"][item["repository_id"]],
                "source_files": item["source_files"],
                "source_lines": item["source_lines"],
                "source_bytes": item["source_bytes"],
                "tokens_M002": item["full_source_tokens"]["M002"],
                "tokens_M003": item["full_source_tokens"]["M003"],
                "M002_context_multiples": item["full_source_tokens"]["M002"] / 65536,
                "M003_context_multiples": item["full_source_tokens"]["M003"] / 65536,
            }
        )
    return rows


def _escape_latex(value: Any) -> str:
    text = str(value)
    for old, new in (
        ("\\", r"\textbackslash{}"),
        ("&", r"\&"),
        ("%", r"\%"),
        ("_", r"\_"),
        ("#", r"\#"),
    ):
        text = text.replace(old, new)
    return text


def paper_tables(
    path: Path,
    repositories: Sequence[dict[str, Any]],
    summaries: Sequence[dict[str, Any]],
    contrast_rows: Sequence[dict[str, Any]],
    reliability: dict[str, Any],
) -> None:
    def reliability_stat(key: str) -> str:
        value = reliability[key]
        return "undefined" if value is None else f"{value:.3f}"

    lines = [
        "% Generated by scripts/analyze_study2.py; do not edit by hand.",
        r"\begin{table}[t]",
        r"\centering\small",
        r"\caption{Study 2 repository scale. Token counts are tokenizer-specific estimates.}",
        r"\label{tab:e08-repositories}",
        r"\begin{tabular}{llrrrr}",
        r"\toprule",
        r"Repo. & Lang. & Tasks & Files & Lines & Qwen tokens \\",
        r"\midrule",
    ]
    for row in repositories:
        lines.append(
            f"{row['repository_id']} & {_escape_latex(row['language'])} & {row['tasks']} & "
            f"{row['source_files']:,} & {row['source_lines']:,} & {row['tokens_M002']:,} \\\\"
        )
    lines.extend([r"\bottomrule", r"\end{tabular}", r"\end{table}", ""])
    lines.extend(
        [
            r"\begin{table*}[t]",
            r"\centering\small",
            r"\caption{Temperature-zero Study 2 outcomes by treatment and model. Intervals are task-bootstrap 95\% intervals.}",
            r"\label{tab:e08-main}",
            r"\begin{tabular}{llrrrrrr}",
            r"\toprule",
            r"Treatment & Model & Resolved & Rate [95\% CI] & Apply & Exact files & Mean tokens & Mean seconds \\",
            r"\midrule",
        ]
    )
    for row in summaries:
        lines.append(
            f"{row['treatment_id']} & {row['model_id']} & {row['resolved_count']}/{row['n']} & "
            f"{row['resolved_rate']:.3f} [{row['resolved_ci_low']:.3f}, {row['resolved_ci_high']:.3f}] & "
            f"{row['patch_apply_rate']:.3f} & {row['exact_modified_file_match_rate']:.3f} & "
            f"{row['mean_total_tokens']:.0f} & {row['mean_elapsed_seconds']:.1f} \\\\"
        )
    lines.extend([r"\bottomrule", r"\end{tabular}", r"\end{table*}", ""])
    lines.extend(
        [
            r"\begin{table*}[t]",
            r"\centering\scriptsize",
            r"\caption{Prespecified paired Study 2 contrasts. Only P1 is confirmatory; the 13 S contrasts form one Holm family.}",
            r"\label{tab:e08-contrasts}",
            r"\begin{tabular}{lllrrrr}",
            r"\toprule",
            r"Contrast & Model & Left--right & RD [95\% CI] & $n_{10}/n_{01}$ & Raw $p$ & Holm $p$ \\",
            r"\midrule",
        ]
    )
    for row in contrast_rows:
        adjusted = "--" if row["family"] == "confirmatory_primary" else f"{row['mcnemar_p_holm']:.4g}"
        lines.append(
            f"{_escape_latex(row['contrast'])} & {row['model_id']} & "
            f"{row['left_treatment']}--{row['right_treatment']} & "
            f"{row['paired_risk_difference']:+.3f} [{row['risk_difference_ci_low']:+.3f}, {row['risk_difference_ci_high']:+.3f}] & "
            f"{row['discordant_left_only']}/{row['discordant_right_only']} & "
            f"{row['mcnemar_p']:.4g} & {adjusted} \\\\"
        )
    lines.extend([r"\bottomrule", r"\end{tabular}", r"\end{table*}", ""])
    lines.extend(
        [
            r"\begin{table}[t]",
            r"\centering\small",
            r"\caption{Temperature-0.2 three-seed reliability sensitivity (24 matched cells; not pooled with the main study).}",
            r"\label{tab:e08-reliability}",
            r"\begin{tabular}{lr}",
            r"\toprule",
            r"Statistic & Value \\",
            r"\midrule",
            f"Unanimous groups & {reliability['unanimous_groups']}/{reliability['groups']} \\\\ ",
            f"Mean pairwise agreement & {reliability['mean_pairwise_agreement']:.3f} \\\\ ",
            f"Fleiss $\\kappa$ & {reliability_stat('fleiss_kappa')} \\\\ ",
            f"Krippendorff $\\alpha$ & {reliability_stat('krippendorff_alpha_nominal')} \\\\ ",
            f"ICC(1,1) & {reliability_stat('icc_1_1')} \\\\ ",
            f"Majority agreement with $T=0$ & {reliability['majority_agreement_with_temperature_zero']:.3f} \\\\ ",
            r"\bottomrule",
            r"\end{tabular}",
            r"\end{table}",
            "",
        ]
    )
    path.write_text("\n".join(lines), encoding="utf-8")


def plots(
    output: Path,
    summaries: Sequence[dict[str, Any]],
    contrast_rows: Sequence[dict[str, Any]],
    strata: Sequence[dict[str, Any]],
    reliability_groups: Sequence[dict[str, Any]],
) -> list[Path]:
    colors = {"M002": "#2364aa", "M003": "#f18f01"}
    x = np.arange(len(TREATMENTS))
    width = 0.36
    fig, axis = plt.subplots(figsize=(8.0, 4.2))
    for offset, model in ((-width / 2, "M002"), (width / 2, "M003")):
        values = [next(r for r in summaries if r["treatment_id"] == t and r["model_id"] == model) for t in TREATMENTS]
        rates = [r["resolved_rate"] for r in values]
        errors = [
            [rate - r["resolved_ci_low"] for rate, r in zip(rates, values)],
            [r["resolved_ci_high"] - rate for rate, r in zip(rates, values)],
        ]
        axis.bar(x + offset, rates, width, label=model, color=colors[model], alpha=0.88)
        axis.errorbar(x + offset, rates, yerr=errors, fmt="none", ecolor="black", capsize=2)
    axis.set_xticks(x, TREATMENTS)
    axis.set(ylabel="Resolved@1", xlabel="Treatment", ylim=(0, 1.0))
    axis.grid(axis="y", alpha=0.25)
    axis.legend(frameon=False)
    fig.tight_layout()
    success_pdf = output / "e08_success_by_treatment_model.pdf"
    fig.savefig(success_pdf)
    fig.savefig(output / "e08_success_by_treatment_model.png", dpi=240)
    plt.close(fig)

    fig, axis = plt.subplots(figsize=(8.0, 6.0))
    forest = list(reversed(contrast_rows))
    y = np.arange(len(forest))
    estimates = np.array([r["paired_risk_difference"] for r in forest])
    low = np.array([r["risk_difference_ci_low"] for r in forest])
    high = np.array([r["risk_difference_ci_high"] for r in 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, [r["contrast"] for r in forest], fontsize=8)
    axis.set(xlabel="Paired resolved@1 risk difference", xlim=(-0.55, 0.55))
    axis.grid(axis="x", alpha=0.25)
    fig.tight_layout()
    forest_pdf = output / "e08_contrast_forest.pdf"
    fig.savefig(forest_pdf)
    fig.savefig(output / "e08_contrast_forest.png", dpi=240)
    plt.close(fig)

    fig, axes = plt.subplots(1, 2, figsize=(10.2, 4.2), sharey=True)
    for axis, model in zip(axes, MODELS):
        matrix = np.array(
            [
                [
                    next(
                        r["resolved_rate"]
                        for r in strata
                        if r["repository_id"] == repo
                        and r["treatment_id"] == treatment
                        and r["model_id"] == model
                    )
                    for treatment in TREATMENTS
                ]
                for repo in REPOSITORIES
            ]
        )
        image = axis.imshow(matrix, vmin=0, vmax=1, cmap="Blues", aspect="auto")
        axis.set_title(model)
        axis.set_xticks(range(len(TREATMENTS)), TREATMENTS, rotation=45, ha="right")
        axis.set_yticks(range(len(REPOSITORIES)), REPOSITORIES)
        for i in range(matrix.shape[0]):
            for j in range(matrix.shape[1]):
                axis.text(j, i, f"{matrix[i, j]:.2f}", ha="center", va="center", fontsize=7,
                          color="white" if matrix[i, j] > 0.55 else "black")
    fig.subplots_adjust(left=0.07, right=0.88, bottom=0.2, top=0.88, wspace=0.10)
    color_axis = fig.add_axes([0.91, 0.20, 0.016, 0.68])
    fig.colorbar(image, cax=color_axis, label="Resolved@1")
    strata_pdf = output / "e08_repository_strata.pdf"
    fig.savefig(strata_pdf)
    fig.savefig(output / "e08_repository_strata.png", dpi=240)
    plt.close(fig)

    distribution = Counter(r["successes_across_seeds"] for r in reliability_groups)
    fig, axis = plt.subplots(figsize=(5.8, 3.8))
    axis.bar(range(4), [distribution.get(i, 0) for i in range(4)], color="#6a4c93")
    axis.set_xticks(range(4))
    axis.set(xlabel="Successes across three temperature-0.2 seeds", ylabel="Matched cells")
    axis.grid(axis="y", alpha=0.25)
    fig.tight_layout()
    reliability_pdf = output / "e08_reliability_distribution.pdf"
    fig.savefig(reliability_pdf)
    fig.savefig(output / "e08_reliability_distribution.png", dpi=240)
    plt.close(fig)
    return [success_pdf, forest_pdf, strata_pdf, reliability_pdf]


def public_cell(row: dict[str, Any]) -> dict[str, Any]:
    return {key: value for key, value in row.items() if key != "protocol_violations"}


def analyze(root: Path) -> dict[str, Any]:
    raw, audit = discover(root)
    main = [row for row in raw if row["repetition"] == 0]
    reliability = [row for row in raw if row["repetition"] == 1]
    summaries = summarize_main(main)
    strata = repository_strata(main)
    contrast_rows = contrasts(main)
    hierarchy = hierarchical_model(main)
    reliability_groups, reliability_summary = reliability_analysis(reliability, main)
    failure_rows, failure_summary = failure_taxonomy(main)
    tool_protocol_rows, tool_protocol_summary = exploratory_tool_protocol_diagnostic(root)
    repositories = repository_characteristics(root)
    output = root / "results/derived/study2"
    output.mkdir(parents=True, exist_ok=True)

    csv_outputs: list[Path] = []
    for filename, rows in (
        ("e08_all_cells.csv", [public_cell(r) for r in raw]),
        ("e08_main_cells.csv", [public_cell(r) for r in main]),
        ("e08_reliability_cells.csv", [public_cell(r) for r in reliability]),
        ("e08_treatment_model_summary.csv", summaries),
        ("e08_repository_strata.csv", strata),
        ("e08_contrasts.csv", contrast_rows),
        ("e08_hierarchical_coefficients.csv", hierarchy["coefficients"]),
        ("e08_reliability_groups.csv", reliability_groups),
        ("e08_failure_taxonomy.csv", failure_rows),
        ("e08_tool_protocol_diagnostic.csv", tool_protocol_rows),
        ("e08_repository_characteristics.csv", repositories),
    ):
        path = output / filename
        write_csv(path, rows)
        csv_outputs.append(path)

    latex_path = output / "e08_paper_tables.tex"
    paper_tables(latex_path, repositories, summaries, contrast_rows, reliability_summary)
    figure_outputs = plots(output, summaries, contrast_rows, strata, reliability_groups)
    analysis_revision = git_output(root, "rev-parse", "HEAD")
    analysis_script = root / "scripts/analyze_study2.py"
    manifest = {
        "schema_version": 1,
        "experiment_id": "E08",
        "execution_revision": EXPECTED_REVISION,
        "analysis_code_revision": analysis_revision,
        "analysis_script_sha256": sha256_file(analysis_script),
        "input_cells": len(raw),
        "input_run_ids": sorted(row["run_id"] for row in raw),
        "raw_manifest_and_metrics_sha256": audit["raw_manifest_and_metrics_sha256"],
        "bootstrap_samples": BOOTSTRAPS,
        "bootstrap_seed": BOOTSTRAP_SEED,
        "primary_test": "two-sided exact McNemar/binomial, alpha=0.05",
        "secondary_multiplicity": "Holm across 13 prespecified model-stratified tests",
        "hierarchical_model": hierarchy["model_type"],
        "reliability_policy": "separate temperature-0.2 sensitivity; never pooled with main",
    }
    manifest_path = output / "analysis_manifest.json"
    manifest_path.write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n", encoding="utf-8")
    report = {
        "schema_version": 1,
        "experiment_id": "E08",
        "claim_boundary": {
            "confirmatory_primary": "H007 versus H000 on M002 only",
            "secondary_family": "13 model-stratified exact paired contrasts with Holm correction",
            "controlled_systems": "A001/A002 are controlled adaptations, not official implementations",
            "language_limit": "Python is represented by one repository and is confounded with repository",
        },
        "audit": audit,
        "repository_characteristics": repositories,
        "treatment_model_summaries": summaries,
        "contrasts": contrast_rows,
        "hierarchical_model": hierarchy,
        "reliability": reliability_summary,
        "failure_analysis": failure_summary,
        "exploratory_tool_protocol_diagnostic": tool_protocol_summary,
        "analysis_manifest": manifest,
    }
    report_path = output / "e08_analysis.json"
    report_path.write_text(json.dumps(report, indent=2, sort_keys=True) + "\n", encoding="utf-8")
    generated = [*csv_outputs, latex_path, *figure_outputs, manifest_path, report_path]
    checksums = {path.name: sha256_file(path) for path in generated}
    checksums_path = output / "SHA256SUMS.json"
    checksums_path.write_text(json.dumps(checksums, indent=2, sort_keys=True) + "\n", encoding="utf-8")
    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])
    args = parser.parse_args()
    try:
        result = analyze(args.root.resolve())
    except (Study2AnalysisError, OSError, ValueError, KeyError, subprocess.CalledProcessError) as exc:
        print(f"STUDY 2 ANALYSIS FAILED: {exc}")
        return 1
    print(json.dumps(result, indent=2, sort_keys=True))
    return 0


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