data: extract the released Figure 5 protocol

1 files changed, 210 insertions, 0 deletions

diff --git a/experiments/extract_dillavou_fig5_protocol.py b/experiments/extract_dillavou_fig5_protocol.py
new file mode 100644
index 0000000..4f0e470
--- /dev/null
+++ b/experiments/extract_dillavou_fig5_protocol.py
@@ -0,0 +1,210 @@
+#!/usr/bin/env python3
+"""Extract the released Figure-5 protocol and endpoints from MATLAB objects."""
+
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+import tempfile
+import zipfile
+
+import numpy as np
+
+try:
+    from matio import load_from_mat
+except ImportError as error:  # pragma: no cover - environment guidance
+    raise SystemExit(
+        "Install the public `mat-io` package to decode MATLAB MCOS objects."
+    ) from error
+
+
+CONDITIONS = {
+    ("taskcycle_L_18", 128.0, 0.0, 25.0): "standard",
+    ("taskcycle_L_19", 32.0, 1.0, 200.0): "overclamp",
+}
+
+
+def scalar(properties: dict, name: str) -> float:
+    return float(np.asarray(properties[name]).reshape(-1)[0])
+
+
+def final_classification_error(confusion: np.ndarray) -> float:
+    final = np.asarray(confusion, dtype=float)[:, :, -1]
+    total = float(np.sum(final))
+    if total <= 0.0:
+        raise ValueError("empty final confusion matrix")
+    return float(1.0 - np.trace(final) / total)
+
+
+def extract_record(path: Path) -> dict | None:
+    properties = load_from_mat(
+        str(path), raw_data=True)["experiment"].properties
+    name = str(np.asarray(properties["Name"]).reshape(-1)[0])
+    if name not in {condition[0] for condition in CONDITIONS}:
+        return None
+    if not all(field in properties for field in ("ETA", "NOR", "ALF")):
+        return None
+    key = (
+        name,
+        scalar(properties, "ETA"),
+        scalar(properties, "NOR"),
+        scalar(properties, "ALF"),
+    )
+    method = CONDITIONS.get(key)
+    if method is None:
+        return None
+    node_multiplier = scalar(properties, "NODEMULT")
+    train = np.asarray(properties["TRAIN"], dtype=float)
+    inputs = train[:2] * node_multiplier
+    center = np.mean(inputs, axis=1)
+    diameter = 2.0 * float(np.mean(np.linalg.norm(
+        inputs - center[:, None], axis=0)))
+    horizontal = np.asarray(
+        properties["HorizontalCapacitors"], dtype=float)
+    vertical = np.asarray(
+        properties["VerticalCapacitors"], dtype=float)
+    gate_multiplier = scalar(properties, "GATEMULT")
+    initial_gates = np.concatenate((
+        horizontal[:, :, 0].reshape(-1),
+        vertical[:, :, 0].reshape(-1),
+    )) * gate_multiplier
+    final_gates = np.concatenate((
+        horizontal[:, :, -1].reshape(-1),
+        vertical[:, :, -1].reshape(-1),
+    )) * gate_multiplier
+    train_mse = np.asarray(properties["TrainMSE"], dtype=float)
+    return {
+        "source_file": path.name,
+        "method": method,
+        "experiment_name": name,
+        "input_diameter_v": diameter,
+        "inputs_v": inputs.tolist(),
+        "classes": np.asarray(
+            properties["TRAINCLASSES"], dtype=int).reshape(-1).tolist(),
+        "source_nodes_zero_indexed": np.asarray(
+            properties["SLOC"], dtype=int).tolist(),
+        "target_nodes_zero_indexed": np.asarray(
+            properties["TLOC"], dtype=int).tolist(),
+        "periodic_axes": np.asarray(
+            properties["ISPERIODIC"], dtype=int).reshape(-1).tolist(),
+        "initial_gates_v": initial_gates.tolist(),
+        "final_gates_v": final_gates.tolist(),
+        "final_classification_error": final_classification_error(
+            np.asarray(properties["TrainConfusion"])),
+        "final_hinge_loss_v2": float(train_mse.reshape(-1)[-1]),
+        "cumulative_learning_time_seconds": float(np.sum(
+            np.asarray(properties["LearnTimes"], dtype=float))) / 1e6,
+        "settings": {
+            "eta_over_129": scalar(properties, "ETA") / 129.0,
+            "alpha_microseconds": scalar(properties, "ALF"),
+            "normalized_alpha": bool(scalar(properties, "NOR")),
+            "hinge_buffer_millivolts": scalar(properties, "BUF"),
+            "one_hot_setting": scalar(properties, "HOT"),
+            "epochs": int(scalar(properties, "EPO")),
+        },
+    }
+
+
+def summarize(records: list[dict]) -> list[dict]:
+    summaries = []
+    methods = sorted({record["method"] for record in records})
+    diameters = sorted({record["input_diameter_v"] for record in records})
+    for method in methods:
+        for diameter in diameters:
+            selected = [
+                record for record in records
+                if record["method"] == method
+                and abs(record["input_diameter_v"] - diameter) < 1e-12
+            ]
+            errors = np.asarray([
+                record["final_classification_error"] for record in selected])
+            hinge = np.asarray([
+                record["final_hinge_loss_v2"] for record in selected])
+            summaries.append({
+                "method": method,
+                "input_diameter_v": diameter,
+                "trials": len(selected),
+                "mean_classification_error": float(np.mean(errors)),
+                "standard_error_classification_error": float(
+                    np.std(errors, ddof=1) / np.sqrt(len(errors))),
+                "mean_hinge_loss_v2": float(np.mean(hinge)),
+                "standard_error_hinge_loss_v2": float(
+                    np.std(hinge, ddof=1) / np.sqrt(len(hinge))),
+            })
+    return summaries
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--artifact-root", type=Path, required=True,
+        help="Extracted maguzj-imperfect-learning-physical-systems source tree")
+    parser.add_argument(
+        "--output", type=Path,
+        default=Path("results/physical_bias/dillavou_fig5_protocol.json"))
+    return parser.parse_args()
+
+
+def main() -> None:
+    args = parse_args()
+    archive = args.artifact_root.resolve() / "big network" / "Experiments.zip"
+    if not archive.exists():
+        raise FileNotFoundError(archive)
+    records = []
+    with tempfile.TemporaryDirectory(prefix="dillavou_fig5_") as directory:
+        directory_path = Path(directory)
+        with zipfile.ZipFile(archive) as handle:
+            names = [
+                name for name in handle.namelist()
+                if name.startswith("Experiments/") and name.endswith(".mat")
+            ]
+            for name in names:
+                path = directory_path / Path(name).name
+                path.write_bytes(handle.read(name))
+                record = extract_record(path)
+                if record is not None:
+                    records.append(record)
+    records.sort(key=lambda record: (
+        record["method"],
+        record["input_diameter_v"],
+        record["classes"],
+    ))
+    if len(records) != 80:
+        raise ValueError(f"expected 80 Figure-5 experiments, found {len(records)}")
+    class_patterns = {
+        tuple(record["classes"]) for record in records
+        if record["method"] == "standard"
+        and record["input_diameter_v"] == min(
+            item["input_diameter_v"] for item in records)
+    }
+    report = {
+        "analysis": "released_dillavou_figure5_protocol",
+        "provenance": {
+            "paper": "Dillavou et al., arXiv:2505.22887v2",
+            "zenodo_record": "15692914",
+            "release": "v1.0.1",
+            "source_revision": "71b8d724afc61d041bcfc1a0b2335dd88b3df62f",
+        },
+        "protocol_checks": {
+            "experiments": len(records),
+            "methods": sorted({record["method"] for record in records}),
+            "input_diameters_v": sorted({
+                record["input_diameter_v"] for record in records}),
+            "label_rotations": len(class_patterns),
+            "trials_per_method_diameter": 8,
+            "grid_shape": [4, 4],
+            "edges": 32,
+        },
+        "summary": summarize(records),
+        "experiments": records,
+    }
+    args.output.parent.mkdir(parents=True, exist_ok=True)
+    args.output.write_text(json.dumps(report, indent=2) + "\n")
+    print(json.dumps(report["protocol_checks"], indent=2))
+    print(json.dumps(report["summary"], indent=2))
+    print(f"wrote {args.output}")
+
+
+if __name__ == "__main__":
+    main()