diff options
Diffstat (limited to 'experiments/analyze_physical_bias_p0.py')
| -rw-r--r-- | experiments/analyze_physical_bias_p0.py | 267 |
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() |