results: reproduce physical structured-bias scaling

1 files changed, 267 insertions, 0 deletions

diff --git a/experiments/analyze_physical_bias_p0.py b/experiments/analyze_physical_bias_p0.py
new file mode 100644
index 0000000..381db13
--- /dev/null
+++ b/experiments/analyze_physical_bias_p0.py
@@ -0,0 +1,267 @@
+#!/usr/bin/env python3
+"""Reproduce physical-bias scaling diagnostics from Dillavou et al. data.
+
+This is a descriptive reanalysis of already published data, not a prospective
+confirmation.  It reads only the small-network CSV files from Zenodo record
+15692914 (release v1.0.1), records their hashes, and measures the low-period
+power laws reported in the paper:
+
+* nonzero combined-error plateau as the task-switching period decreases;
+* squared cycle span proportional to approximately period squared;
+* therefore a nonzero cycle-speed proxy span / period.
+
+The script also audits the released two-dimensional bow-tie trajectories.
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+from pathlib import Path
+import re
+from typing import Dict, List
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+ZENODO_RECORD = "15692914"
+ZENODO_DOI = "10.5281/zenodo.15692914"
+RELEASE = "v1.0.1"
+SOURCE_TREE = "maguzj-imperfect-learning-physical-systems-71b8d72"
+
+
+def sha256(path: Path) -> str:
+    digest = hashlib.sha256()
+    with path.open("rb") as handle:
+        for block in iter(lambda: handle.read(1024 * 1024), b""):
+            digest.update(block)
+    return digest.hexdigest()
+
+
+def log_slope(values: np.ndarray) -> float:
+    if values.ndim != 2 or values.shape[1] != 2:
+        raise ValueError("expected a two-column positive-valued table")
+    if np.any(values <= 0):
+        raise ValueError("log slope requires strictly positive values")
+    return float(np.polyfit(np.log(values[:, 0]), np.log(values[:, 1]), 1)[0])
+
+
+def load_table(path: Path) -> np.ndarray:
+    values = np.loadtxt(path, delimiter=",")
+    if values.ndim != 2 or values.shape[1] != 2:
+        raise ValueError(f"unexpected table shape in {path}: {values.shape}")
+    if not np.all(np.isfinite(values)):
+        raise ValueError(f"nonfinite value in {path}")
+    if np.any(np.diff(values[:, 0]) <= 0):
+        raise ValueError(f"periods are not strictly increasing in {path}")
+    return values
+
+
+def experiment_tables(data_dir: Path, fit_points: int) -> Dict[str, dict]:
+    records: Dict[str, dict] = {}
+    for experiment in (1, 2, 3):
+        mse_path = data_dir / f"measured-MSE-exp{experiment}.csv"
+        span2_path = data_dir / f"measured-DG2-exp{experiment}.csv"
+        mse = load_table(mse_path)
+        span2 = load_table(span2_path)
+        if not np.array_equal(mse[:, 0], span2[:, 0]):
+            raise ValueError(f"period grids disagree for experiment {experiment}")
+        if len(mse) < fit_points:
+            raise ValueError("fit_points exceeds a source table")
+        low_mse = mse[:fit_points]
+        low_span2 = span2[:fit_points]
+        speed = np.sqrt(low_span2[:, 1]) / low_span2[:, 0]
+        records[str(experiment)] = {
+            "period": mse[:, 0].tolist(),
+            "combined_error": mse[:, 1].tolist(),
+            "cycle_span_squared": span2[:, 1].tolist(),
+            "low_period_points": int(fit_points),
+            "low_period_combined_error_log_slope": log_slope(low_mse),
+            "low_period_span_squared_log_slope": log_slope(low_span2),
+            "low_period_span_log_slope": 0.5 * log_slope(low_span2),
+            "low_period_speed_proxy": speed.tolist(),
+            "low_period_speed_proxy_mean": float(np.mean(speed)),
+            "low_period_speed_proxy_cv": float(np.std(speed) / np.mean(speed)),
+            "source_files": {
+                "combined_error": {
+                    "path": str(mse_path),
+                    "sha256": sha256(mse_path),
+                },
+                "cycle_span_squared": {
+                    "path": str(span2_path),
+                    "sha256": sha256(span2_path),
+                },
+            },
+        }
+    return records
+
+
+def bowtie_records(bowtie_dir: Path) -> List[dict]:
+    pattern = re.compile(r"bowtie-exp(?P<experiment>[123])-(?P<period>[0-9.]+)s\.csv")
+    records = []
+    for path in sorted(bowtie_dir.glob("bowtie-exp*-*s.csv")):
+        match = pattern.fullmatch(path.name)
+        if match is None:
+            raise ValueError(f"unrecognized bow-tie filename: {path.name}")
+        trajectory = np.loadtxt(path, delimiter=",")
+        if trajectory.ndim != 2 or trajectory.shape[0] != 2:
+            raise ValueError(f"unexpected bow-tie shape in {path}: {trajectory.shape}")
+        pairwise = trajectory[:, :, None] - trajectory[:, None, :]
+        diameter = float(np.sqrt(np.sum(pairwise * pairwise, axis=0)).max())
+        midpoint = trajectory.shape[1] // 2
+        half_cycle_span = float(np.linalg.norm(
+            trajectory[:, 0] - trajectory[:, midpoint]))
+        path_length = float(np.linalg.norm(
+            np.diff(trajectory, axis=1), axis=0).sum())
+        records.append({
+            "experiment": int(match.group("experiment")),
+            "period_seconds": float(match.group("period")),
+            "samples": int(trajectory.shape[1]),
+            "diameter": diameter,
+            "half_cycle_span": half_cycle_span,
+            "path_length": path_length,
+            "trajectory": trajectory.tolist(),
+            "source_file": {"path": str(path), "sha256": sha256(path)},
+        })
+    if len(records) != 10:
+        raise ValueError(f"expected 10 bow-tie trajectories, found {len(records)}")
+    return records
+
+
+def representative_quadratic(period: np.ndarray, values: np.ndarray) -> np.ndarray:
+    return values[0] * (period / period[0]) ** 2
+
+
+def plot_report(report: dict, output: Path) -> None:
+    colors = ["#4477AA", "#EE6677", "#228833"]
+    fig, axes = plt.subplots(2, 2, figsize=(9.0, 7.0))
+
+    for index, (experiment, record) in enumerate(report["experiments"].items()):
+        period = np.asarray(record["period"])
+        error = np.asarray(record["combined_error"])
+        span2 = np.asarray(record["cycle_span_squared"])
+        label = f"physical experiment {experiment}"
+        axes[0, 0].loglog(period, error, "o-", color=colors[index], label=label)
+        axes[0, 1].loglog(period, span2, "o-", color=colors[index], label=label)
+        low_n = int(record["low_period_points"])
+        axes[0, 1].loglog(
+            period[:low_n], representative_quadratic(period[:low_n], span2[:low_n]),
+            "--", color=colors[index], alpha=0.55)
+        speed = np.sqrt(span2) / period
+        axes[1, 0].semilogx(period, speed, "o-", color=colors[index], label=label)
+
+    axes[0, 0].set_title("A  Error remains nonzero under rapid switching")
+    axes[0, 0].set_xlabel("task-switching period (s)")
+    axes[0, 0].set_ylabel("combined error")
+    axes[0, 0].legend(frameon=False, fontsize=8)
+
+    axes[0, 1].set_title("B  Cycle span follows the bias-drift prediction")
+    axes[0, 1].set_xlabel("task-switching period (s)")
+    axes[0, 1].set_ylabel("cycle span squared")
+
+    axes[1, 0].set_title("C  Span per unit period does not vanish")
+    axes[1, 0].set_xlabel("task-switching period (s)")
+    axes[1, 0].set_ylabel("sqrt(span squared) / period")
+
+    selected = [
+        item for item in report["bowtie_trajectories"]
+        if item["experiment"] == 1
+    ]
+    selected.sort(key=lambda item: item["period_seconds"])
+    for index, item in enumerate(selected):
+        trajectory = np.asarray(item["trajectory"])
+        axes[1, 1].plot(
+            trajectory[0], trajectory[1], "o-", markersize=2.1,
+            linewidth=1.0, color=colors[index],
+            label=f"period={item['period_seconds']:g} s")
+    axes[1, 1].set_title("D  Measured bias-driven parameter cycles")
+    axes[1, 1].set_xlabel("gate voltage +")
+    axes[1, 1].set_ylabel("gate voltage -")
+    axes[1, 1].legend(frameon=False, fontsize=8)
+
+    for axis in axes.flat:
+        axis.grid(alpha=0.18)
+    fig.suptitle(
+        "Published physical coupled-learning data (Dillavou et al.; Zenodo 15692914)",
+        fontsize=11)
+    fig.tight_layout()
+    output.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(output, dpi=180)
+    plt.close(fig)
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--artifact-root", type=Path, required=True,
+        help="Root of the extracted Zenodo source tree")
+    parser.add_argument("--fit-points", type=int, default=6)
+    parser.add_argument(
+        "--json", type=Path,
+        default=Path("results/physical_bias/p0_summary.json"))
+    parser.add_argument(
+        "--figure", type=Path,
+        default=Path("results/figs/physical_bias_p0.png"))
+    return parser.parse_args()
+
+
+def main() -> None:
+    args = parse_args()
+    if args.fit_points < 3:
+        raise ValueError("fit_points must be at least three")
+    root = args.artifact_root.resolve()
+    expected_name = SOURCE_TREE
+    if root.name != expected_name:
+        raise ValueError(f"expected source root {expected_name}, received {root.name}")
+    small = root / "small network"
+    records = experiment_tables(small / "MSE-DG2-Data", args.fit_points)
+    bowties = bowtie_records(small / "experimental_bowties")
+    mse_slopes = np.asarray([
+        value["low_period_combined_error_log_slope"]
+        for value in records.values()])
+    span2_slopes = np.asarray([
+        value["low_period_span_squared_log_slope"]
+        for value in records.values()])
+    report = {
+        "analysis": "published_physical_bias_descriptive_reproduction",
+        "confirmatory": False,
+        "provenance": {
+            "zenodo_record": ZENODO_RECORD,
+            "doi": ZENODO_DOI,
+            "release": RELEASE,
+            "source_tree": SOURCE_TREE,
+        },
+        "fit_definition": {
+            "points": int(args.fit_points),
+            "selection": "lowest task-switching periods in each published table",
+            "regression": "ordinary least squares on log(period), log(metric)",
+        },
+        "experiments": records,
+        "bowtie_trajectories": bowties,
+        "summary": {
+            "combined_error_log_slope_mean": float(np.mean(mse_slopes)),
+            "combined_error_log_slopes": mse_slopes.tolist(),
+            "span_squared_log_slope_mean": float(np.mean(span2_slopes)),
+            "span_squared_log_slopes": span2_slopes.tolist(),
+            "span_log_slope_mean": float(0.5 * np.mean(span2_slopes)),
+            "descriptive_plateau_check_abs_error_slope_below_0p4": bool(
+                np.all(np.abs(mse_slopes) < 0.4)),
+            "descriptive_linear_drift_check_span2_slope_1p5_to_2p3": bool(
+                np.all((span2_slopes > 1.5) & (span2_slopes < 2.3))),
+        },
+    }
+    args.json.parent.mkdir(parents=True, exist_ok=True)
+    args.json.write_text(json.dumps(report, indent=2) + "\n")
+    plot_report(report, args.figure)
+    print(json.dumps(report["summary"], indent=2))
+    print(f"wrote {args.json}")
+    print(f"wrote {args.figure}")
+
+
+if __name__ == "__main__":
+    main()