#!/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())