analysis: freeze constant-calibration stress comparison

1 files changed, 222 insertions, 0 deletions

diff --git a/experiments/analyze_constant_calibration_stress.py b/experiments/analyze_constant_calibration_stress.py
new file mode 100644
index 0000000..4eeb405
--- /dev/null
+++ b/experiments/analyze_constant_calibration_stress.py
@@ -0,0 +1,222 @@
+#!/usr/bin/env python3
+"""Analyze the frozen side-32 static-calibration budget stress test."""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+from pathlib import Path
+
+import numpy as np
+
+
+TASKS = tuple(range(40))
+SEEDS = (20260830, 20260831, 20260832)
+SIDE = 32
+EDGES = 2048
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--core",
+        type=Path,
+        default=Path("results/coupled_ladder/p2_confirm_side32.json"),
+    )
+    parser.add_argument(
+        "--stress",
+        type=Path,
+        default=Path("results/coupled_ladder/p4_constant256_side32.json"),
+    )
+    parser.add_argument(
+        "--output",
+        type=Path,
+        default=Path("results/coupled_ladder/p4_constant256_analysis.json"),
+    )
+    parser.add_argument(
+        "--output-csv",
+        type=Path,
+        default=Path("results/coupled_ladder/p4_constant256_source.csv"),
+    )
+    parser.add_argument("--bootstrap-replicates", type=int, default=20000)
+    parser.add_argument("--bootstrap-seed", type=int, default=20260829)
+    return parser.parse_args()
+
+
+def validate(report: dict, observations: int, methods: set[str]) -> None:
+    protocol = report["protocol"]
+    if not report.get("confirmatory", False):
+        raise ValueError("source must be labeled confirmatory")
+    if protocol["sizes"] != [SIDE]:
+        raise ValueError("source must contain only side 32")
+    if protocol["rotations_per_input_diameter"] != 8:
+        raise ValueError("source must contain all eight rotations")
+    if tuple(protocol["device_seeds"]) != SEEDS:
+        raise ValueError("source device seeds do not match")
+    if set(protocol["methods"]) != methods:
+        raise ValueError("source methods do not match")
+    if protocol["epochs"] != 600:
+        raise ValueError("source must use 600 epochs")
+    if abs(float(protocol["learning_time_seconds_by_side"]["32"]) - 0.01) > 1e-15:
+        raise ValueError("source learning exposure does not match")
+    if protocol["calibration_observations"] != observations:
+        raise ValueError("source calibration budget does not match")
+    records = report["records"]
+    expected = {(SIDE, task, seed) for task in TASKS for seed in SEEDS}
+    actual = {
+        (record["side"], record["task_index"], record["device_seed"])
+        for record in records
+    }
+    if actual != expected:
+        raise ValueError("source does not contain the exact 120 paired cells")
+    if any(set(record["methods"]) != methods for record in records):
+        raise ValueError("each source cell must contain exactly its protocol methods")
+    if any(
+        value["status"] != "completed"
+        for record in records
+        for value in record["methods"].values()
+    ):
+        raise ValueError("all source methods must complete")
+
+
+def task_values(report: dict, method: str, metric: str) -> np.ndarray:
+    values = []
+    for task in TASKS:
+        device_values = [
+            record["methods"][method][metric]
+            for record in report["records"]
+            if record["task_index"] == task
+        ]
+        values.append(float(np.mean(device_values)))
+    return np.asarray(values)
+
+
+def paired_summary(
+    budget16: np.ndarray,
+    budget256: np.ndarray,
+    *,
+    higher_is_better: bool,
+    rng: np.random.Generator,
+    replicates: int,
+) -> dict:
+    if higher_is_better:
+        improvement = budget256 - budget16
+    else:
+        improvement = budget16 - budget256
+    indices = rng.integers(0, len(improvement), size=(replicates, len(improvement)))
+    draws = np.mean(improvement[indices], axis=1)
+    return {
+        "budget16_mean": float(np.mean(budget16)),
+        "budget256_mean": float(np.mean(budget256)),
+        "budget256_improvement": float(np.mean(improvement)),
+        "paired_task_bootstrap_95ci": [
+            float(value) for value in np.percentile(draws, (2.5, 97.5))
+        ],
+        "higher_is_better": higher_is_better,
+    }
+
+
+def observation_cost(report: dict, observations: int) -> np.ndarray:
+    output = []
+    for task in TASKS:
+        device_values = []
+        for record in report["records"]:
+            if record["task_index"] != task:
+                continue
+            restricted_edges = record["methods"]["constant"][
+                "restricted_edge_updates_to_stable_zero_error"
+            ]
+            device_values.append(restricted_edges + observations * EDGES)
+        output.append(float(np.mean(device_values)))
+    return np.asarray(output)
+
+
+def main() -> None:
+    args = parse_args()
+    core = json.loads(args.core.read_text())
+    stress = json.loads(args.stress.read_text())
+    validate(
+        core,
+        observations=16,
+        methods={"clean", "matched_noise", "raw", "constant", "sdil"},
+    )
+    validate(stress, observations=256, methods={"constant"})
+    rng = np.random.default_rng(args.bootstrap_seed)
+    metrics = {
+        "classification_error": paired_summary(
+            task_values(core, "constant", "classification_error"),
+            task_values(stress, "constant", "classification_error"),
+            higher_is_better=False,
+            rng=rng,
+            replicates=args.bootstrap_replicates,
+        ),
+        "classification_error_auc": paired_summary(
+            task_values(core, "constant", "classification_error_auc"),
+            task_values(stress, "constant", "classification_error_auc"),
+            higher_is_better=False,
+            rng=rng,
+            replicates=args.bootstrap_replicates,
+        ),
+        "stable_zero_fraction": paired_summary(
+            task_values(core, "constant", "reached_stable_zero_error"),
+            task_values(stress, "constant", "reached_stable_zero_error"),
+            higher_is_better=True,
+            rng=rng,
+            replicates=args.bootstrap_replicates,
+        ),
+        "local_scalar_reads_to_target": paired_summary(
+            observation_cost(core, 16),
+            observation_cost(stress, 256),
+            higher_is_better=False,
+            rng=rng,
+            replicates=args.bootstrap_replicates,
+        ),
+    }
+    error_result = metrics["classification_error"]
+    sampling_limited = bool(
+        error_result["budget256_improvement"] >= 0.02
+        and error_result["paired_task_bootstrap_95ci"][0] > 0.0
+    )
+    analysis = {
+        "analysis": "side32_static_calibration_budget_stress",
+        "confirmatory": True,
+        "comparison": "16 versus 256 instruction-off observations per edge",
+        "bootstrap": {
+            "unit": "task; three component draws averaged within task",
+            "task_clusters": 40,
+            "replicates": args.bootstrap_replicates,
+            "seed": args.bootstrap_seed,
+            "interval": "percentile 95%",
+        },
+        "metrics": metrics,
+        "decision_rule": (
+            "sampling-limited only if final error improves by at least 0.02 "
+            "and the paired 95% interval excludes zero"
+        ),
+        "sampling_limited": sampling_limited,
+        "sources": {"budget16": str(args.core), "budget256": str(args.stress)},
+    }
+    args.output.parent.mkdir(parents=True, exist_ok=True)
+    args.output.write_text(json.dumps(analysis, indent=2) + "\n")
+    args.output_csv.parent.mkdir(parents=True, exist_ok=True)
+    with args.output_csv.open("w", newline="") as stream:
+        writer = csv.writer(stream)
+        writer.writerow((
+            "metric", "budget16_mean", "budget256_mean", "improvement",
+            "ci_low", "ci_high",
+        ))
+        for metric, values in metrics.items():
+            writer.writerow((
+                metric,
+                values["budget16_mean"],
+                values["budget256_mean"],
+                values["budget256_improvement"],
+                values["paired_task_bootstrap_95ci"][0],
+                values["paired_task_bootstrap_95ci"][1],
+            ))
+    print(json.dumps(analysis, indent=2))
+
+
+if __name__ == "__main__":
+    main()