diff --git a/examples/paper_runs/paper_04/03_analyze_comparison.sh b/examples/paper_runs/paper_04/03_analyze_comparison.sh index 551f8ff..f60e3ba 100755 --- a/examples/paper_runs/paper_04/03_analyze_comparison.sh +++ b/examples/paper_runs/paper_04/03_analyze_comparison.sh @@ -38,4 +38,12 @@ done < <(paper_resolve_decoders) --out-md "${OUT_DIR}/table_source_vs_lidmas.md" \ --out-prefix "${OUT_DIR}/figure_source_vs_lidmas" +"${PY_BIN}" "${SCRIPT_DIR}/scripts/analyze_logical_error_rate.py" \ + --truth-dir "${REQUEST_DIR}" \ + --replay-manifest "${REPLAY_DIR}/replay_manifest.csv" \ + --responses-dir "${REPLAY_DIR}" \ + --out-csv "${OUT_DIR}/table_logical_error_rate.csv" \ + --out-md "${OUT_DIR}/table_logical_error_rate.md" \ + --out-prefix "${OUT_DIR}/figure_logical_error_rate" + echo "paper_04 step 03 complete: ${OUT_DIR}" diff --git a/examples/paper_runs/paper_04/08_code_family_comparison.sh b/examples/paper_runs/paper_04/08_code_family_comparison.sh index 72597f1..1d9b0a6 100755 --- a/examples/paper_runs/paper_04/08_code_family_comparison.sh +++ b/examples/paper_runs/paper_04/08_code_family_comparison.sh @@ -79,4 +79,10 @@ done --manifest "${MANIFEST}" \ --out-dir "${OUT_DIR}" +"${PY_BIN}" "${SCRIPT_DIR}/scripts/compose_journal_results_figure.py" \ + --analysis-dir "${OUT_DIR}" \ + --out-prefix "${OUT_DIR}/figure_journal_results_summary" \ + --manuscript-dir "${OUT_DIR}/manuscript_figures" \ + --write-standalone + echo "paper_04 unified analysis complete: ${OUT_DIR}" diff --git a/examples/paper_runs/paper_04/scripts/analyze_code_family_trends.py b/examples/paper_runs/paper_04/scripts/analyze_code_family_trends.py index 0f0fa3e..2469bef 100755 --- a/examples/paper_runs/paper_04/scripts/analyze_code_family_trends.py +++ b/examples/paper_runs/paper_04/scripts/analyze_code_family_trends.py @@ -16,6 +16,8 @@ def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--manifest", required=True, help="Family run manifest CSV.") parser.add_argument("--out-dir", required=True, help="Output directory.") + parser.add_argument("--rank-bootstrap", type=int, default=4000, help="Bootstrap samples for rank-stability analysis.") + parser.add_argument("--rank-seed", type=int, default=20260410, help="RNG seed for rank-stability bootstrap.") return parser.parse_args() @@ -55,6 +57,10 @@ def _write_csv(rows: list[dict[str, Any]], fields: list[str], out_csv: Path) -> or key.endswith("_low") or key.endswith("_high") or key.endswith("_ratio") + or key.endswith("_share") + or key.endswith("_squares") + or key.endswith("_square") + or key.endswith("_prob") ): if key in out: out[key] = _fmt(out[key]) @@ -63,7 +69,11 @@ def _write_csv(rows: list[dict[str, Any]], fields: list[str], out_csv: Path) -> def _plot_family_tradeoff(summary_rows: list[dict[str, Any]], families: list[str], out_prefix: Path) -> None: try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) import matplotlib.pyplot as plt # type: ignore + from matplotlib.lines import Line2D # type: ignore except Exception as exc: print(f"Warning: matplotlib unavailable; skipping family tradeoff plot ({exc}).") return @@ -72,30 +82,74 @@ def _plot_family_tradeoff(summary_rows: list[dict[str, Any]], families: list[str return cols = len(families) - fig, axes = plt.subplots(1, cols, figsize=(5.0 * max(1, cols), 4.4), dpi=320, constrained_layout=True) + fig, axes = plt.subplots(1, cols, figsize=(5.6 * max(1, cols), 4.3), dpi=320, constrained_layout=False) axes_arr = np.atleast_1d(axes) - palette = ["#1f77b4", "#2ca02c", "#d62728", "#9467bd", "#8c564b", "#17becf"] + palette = {"bp": "#1f77b4", "mwpm": "#2ca02c", "uf": "#d62728"} for i, fam in enumerate(families): ax = axes_arr[i] fam_rows = [r for r in summary_rows if str(r.get("family", "")) == fam] - fam_rows = sorted(fam_rows, key=lambda r: str(r.get("decoder", ""))) + fam_rows = sorted(fam_rows, key=lambda r: _f(r.get("mean_avg_flip_sources"))) decoders = [str(r.get("decoder", "")) for r in fam_rows] x = np.asarray([_f(r.get("mean_avg_flip_sources")) for r in fam_rows], dtype=float) - y = np.asarray([_f(r.get("mean_warning_rate_sources")) for r in fam_rows], dtype=float) + y = np.asarray([_f(r.get("mean_nonempty_flip_rate_sources")) for r in fam_rows], dtype=float) + delta = np.asarray([_f(r.get("mean_abs_source_delta_flip")) for r in fam_rows], dtype=float) + finite_delta = delta[np.isfinite(delta)] + if finite_delta.size: + d_min = float(np.min(finite_delta)) + d_span = float(np.max(finite_delta) - d_min) + sizes = 150.0 + 420.0 * ((delta - d_min) / d_span if d_span > 1e-12 else np.ones_like(delta) * 0.45) + else: + sizes = np.ones_like(x) * 260.0 for j, dec in enumerate(decoders): - ax.scatter(x[j], y[j], s=56, color=palette[j % len(palette)], label=dec) - ax.text(x[j], y[j], dec, fontsize=8, ha="left", va="bottom") + color = palette.get(dec, "#6B7280") + ax.scatter(x[j], y[j], s=sizes[j], color=color, edgecolors="white", linewidths=0.9, zorder=3) + ax.annotate( + dec.upper(), + xy=(x[j], y[j]), + xytext=(5, 8), + textcoords="offset points", + fontsize=8.2, + fontweight="bold", + ha="left", + va="bottom", + color="#111827", + ) - ax.set_title(f"{fam}: Decoder Tradeoff") + if np.isfinite(x).any(): + x_lo = float(np.nanmin(x)) + x_hi = float(np.nanmax(x)) + pad = max(0.08 * (x_hi - x_lo), 0.08) + ax.set_xlim(x_lo - pad, x_hi + pad) + if np.isfinite(y).any(): + y_lo = float(np.nanmin(y)) + y_hi = float(np.nanmax(y)) + pad = max(0.16 * (y_hi - y_lo), 0.015) + ax.set_ylim(max(0.0, y_lo - pad), min(1.04, y_hi + pad)) + + panel = chr(ord("a") + i) + ax.text(-0.12, 1.06, panel, transform=ax.transAxes, fontsize=13, fontweight="bold", va="top") + ax.set_title(f"{fam.upper()} decoder operating point", fontsize=11.5, pad=10) ax.set_xlabel("Mean flip count (sources)") - ax.set_ylabel("Warning rate (sources)") - ax.grid(alpha=0.25) + ax.set_ylabel("Nonempty flip rate (sources)") + ax.grid(alpha=0.25, linewidth=0.8) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) - handles, labels = axes_arr[0].get_legend_handles_labels() - if handles: - fig.legend(handles, labels, loc="upper center", ncol=min(6, len(labels)), frameon=False, bbox_to_anchor=(0.5, 1.05)) + size_handles = [ + Line2D([0], [0], marker="o", color="none", markerfacecolor="#9CA3AF", markeredgecolor="white", markersize=7, label="lower source sensitivity"), + Line2D([0], [0], marker="o", color="none", markerfacecolor="#9CA3AF", markeredgecolor="white", markersize=12, label="higher source sensitivity"), + ] + fig.legend( + handles=size_handles, + loc="lower center", + ncol=2, + frameon=False, + bbox_to_anchor=(0.5, -0.02), + fontsize=8.0, + ) + fig.subplots_adjust(left=0.08, right=0.99, top=0.84, bottom=0.22, wspace=0.28) for ext in (".png", ".pdf", ".svg"): fig.savefig(out_prefix.with_suffix(ext), bbox_inches="tight") @@ -104,6 +158,9 @@ def _plot_family_tradeoff(summary_rows: list[dict[str, Any]], families: list[str def _plot_family_delta_forest(delta_rows: list[dict[str, Any]], families: list[str], out_prefix: Path) -> None: try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) import matplotlib.pyplot as plt # type: ignore except Exception as exc: print(f"Warning: matplotlib unavailable; skipping family delta forest ({exc}).") @@ -148,6 +205,142 @@ def _plot_family_delta_forest(delta_rows: list[dict[str, Any]], families: list[s plt.close(fig) +def _plot_source_vs_lidmas( + source_rows: list[dict[str, Any]], + families: list[str], + out_prefix: Path, +) -> None: + try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt # type: ignore + from matplotlib.lines import Line2D # type: ignore + except Exception as exc: + print(f"Warning: matplotlib unavailable; skipping source-vs-lidmas plot ({exc}).") + return + + if not source_rows: + return + + clean_rows = [] + for row in source_rows: + if str(row.get("status", "")).strip() != "ok": + continue + s = _f(row.get("avg_flip_count_source")) + r = _f(row.get("avg_flip_count_reference")) + if not (np.isfinite(s) and np.isfinite(r)): + continue + clean_rows.append(row) + if not clean_rows: + return + + decoders = sorted({str(r.get("decoder", "")) for r in clean_rows if str(r.get("decoder", ""))}) + sources = sorted({str(r.get("source_dataset", "")) for r in clean_rows if str(r.get("source_dataset", ""))}) + if not decoders or not sources: + return + + src_colors = { + "pennylane": "#1f77b4", + "qiskit": "#2ca02c", + "cirq": "#9467bd", + } + source_offsets = np.linspace(-0.22, 0.22, num=len(sources)) if len(sources) > 1 else np.asarray([0.0]) + + cols = max(1, len(families)) + fig, axes = plt.subplots(1, cols, figsize=(6.2 * cols, 4.8), dpi=320, constrained_layout=True, sharey=True) + axes_arr = np.atleast_1d(axes) + + for i, fam in enumerate(families): + ax = axes_arr[i] + fam_rows = [r for r in clean_rows if str(r.get("family", "")) == fam] + if not fam_rows: + ax.set_axis_off() + continue + + x_centers = np.arange(len(decoders), dtype=float) + for d_idx, dec in enumerate(decoders): + d_rows = [r for r in fam_rows if str(r.get("decoder", "")) == dec] + if not d_rows: + continue + + ref_vals = [_f(r.get("avg_flip_count_reference")) for r in d_rows if np.isfinite(_f(r.get("avg_flip_count_reference")))] + if ref_vals: + ref = float(np.mean(ref_vals)) + ax.scatter( + x_centers[d_idx], + ref, + marker="D", + s=44, + color="#111827", + edgecolors="white", + linewidths=0.6, + zorder=4, + ) + + for s_idx, src in enumerate(sources): + row = next((r for r in d_rows if str(r.get("source_dataset", "")) == src), None) + if row is None: + continue + y_src = _f(row.get("avg_flip_count_source")) + y_ref = _f(row.get("avg_flip_count_reference")) + if not (np.isfinite(y_src) and np.isfinite(y_ref)): + continue + x_src = x_centers[d_idx] + float(source_offsets[s_idx]) + ax.plot([x_centers[d_idx], x_src], [y_ref, y_src], color="#9CA3AF", linewidth=0.9, alpha=0.7, zorder=1) + ax.scatter( + x_src, + y_src, + marker="o", + s=42, + color=src_colors.get(src, "#6B7280"), + edgecolors="white", + linewidths=0.6, + zorder=3, + ) + + ax.set_xticks(x_centers) + ax.set_xticklabels(decoders) + ax.set_xlabel("Decoder") + if i == 0: + ax.set_ylabel("Average flip count") + ax.set_title(f"{fam}: Source vs LiDMaS+ reference") + ax.grid(axis="y", alpha=0.25) + + legend_items: list[Line2D] = [ + Line2D( + [0], + [0], + marker="D", + color="none", + markerfacecolor="#111827", + markeredgecolor="white", + markeredgewidth=0.6, + markersize=7, + label="LiDMaS+ reference", + ) + ] + for src in sources: + legend_items.append( + Line2D( + [0], + [0], + marker="o", + color="none", + markerfacecolor=src_colors.get(src, "#6B7280"), + markeredgecolor="white", + markeredgewidth=0.6, + markersize=7, + label=src, + ) + ) + fig.legend(handles=legend_items, loc="upper center", ncol=min(5, len(legend_items)), frameon=False, bbox_to_anchor=(0.5, 1.06)) + + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_prefix.with_suffix(ext), bbox_inches="tight") + plt.close(fig) + + def _plot_normalized_trends( norm_rows: list[dict[str, Any]], families: list[str], @@ -155,6 +348,9 @@ def _plot_normalized_trends( out_prefix: Path, ) -> None: try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) import matplotlib.pyplot as plt # type: ignore except Exception as exc: print(f"Warning: matplotlib unavailable; skipping normalized trends plot ({exc}).") @@ -201,6 +397,392 @@ def _plot_normalized_trends( plt.close(fig) +def _rank_stability_by_family( + matrix_rows: list[dict[str, Any]], + families: list[str], + decoders: list[str], + *, + bootstrap: int, + seed: int, +) -> tuple[list[dict[str, Any]], dict[str, np.ndarray], list[int]]: + ok_rows = [r for r in matrix_rows if str(r.get("status", "")).strip() == "ok"] + ranks = list(range(1, len(decoders) + 1)) + rng = np.random.default_rng(seed) + out_rows: list[dict[str, Any]] = [] + heatmaps: dict[str, np.ndarray] = {} + + for fam in families: + fam_rows = [r for r in ok_rows if str(r.get("family", "")) == fam] + sources = sorted( + { + str(r.get("dataset", "")).strip() + for r in fam_rows + if str(r.get("dataset", "")).strip() and str(r.get("dataset", "")).strip() != "lidmas_reference" + } + ) + probs = np.zeros((len(decoders), len(decoders)), dtype=float) + if not sources or not decoders: + heatmaps[fam] = probs + continue + + metric_map: dict[tuple[str, str], float] = {} + for row in fam_rows: + source = str(row.get("dataset", "")).strip() + decoder = str(row.get("decoder", "")).strip() + if source and decoder: + metric_map[(source, decoder)] = _f(row.get("avg_flip_count")) + + counts = np.zeros_like(probs) + source_idx = np.arange(len(sources), dtype=int) + n_boot = max(1, int(bootstrap)) + for _ in range(n_boot): + sample = rng.choice(source_idx, size=len(source_idx), replace=True) + means: list[float] = [] + for decoder in decoders: + vals = [ + metric_map.get((sources[int(idx)], decoder), float("nan")) + for idx in sample + ] + finite = [v for v in vals if np.isfinite(v)] + means.append(float(np.mean(finite)) if finite else float("inf")) + order = sorted(range(len(decoders)), key=lambda idx: (means[idx], decoders[idx])) + for rank_zero, decoder_idx in enumerate(order): + counts[decoder_idx, rank_zero] += 1.0 + + probs = counts / float(n_boot) + heatmaps[fam] = probs + for decoder_idx, decoder in enumerate(decoders): + for rank_zero, rank in enumerate(ranks): + out_rows.append( + { + "family": fam, + "decoder": decoder, + "rank": rank, + "rank_prob": float(probs[decoder_idx, rank_zero]), + "bootstrap_samples": n_boot, + "source_count": len(sources), + } + ) + + return out_rows, heatmaps, ranks + + +def _plot_rank_stability_by_family( + heatmaps: dict[str, np.ndarray], + families: list[str], + decoders: list[str], + ranks: list[int], + out_prefix: Path, +) -> None: + try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt # type: ignore + except Exception as exc: + print(f"Warning: matplotlib unavailable; skipping rank-stability plot ({exc}).") + return + + if not heatmaps or not decoders or not ranks: + return + + cols = max(1, len(families)) + fig, axes = plt.subplots(1, cols, figsize=(5.4 * cols, 3.6), dpi=320, constrained_layout=False) + axes_arr = np.atleast_1d(axes) + + last_im = None + for i, fam in enumerate(families): + ax = axes_arr[i] + probs = heatmaps.get(fam, np.zeros((len(decoders), len(ranks)), dtype=float)) + last_im = ax.imshow(probs, cmap="YlGnBu", vmin=0.0, vmax=1.0, aspect="auto") + ax.set_title(f"{fam.upper()} rank stability", fontsize=11.5, pad=9) + ax.set_xticks(np.arange(len(ranks))) + ax.set_xticklabels([str(rank) for rank in ranks]) + ax.set_xlabel("Rank (1 = lowest flip count)") + ax.set_yticks(np.arange(len(decoders))) + ax.set_yticklabels([decoder.upper() for decoder in decoders]) + ax.text(-0.14, 1.08, chr(ord("a") + i), transform=ax.transAxes, fontsize=13, fontweight="bold", va="top") + + for row in range(probs.shape[0]): + for col in range(probs.shape[1]): + value = float(probs[row, col]) + color = "white" if value > 0.62 else "#111827" + ax.text(col, row, f"{value:.2f}", ha="center", va="center", fontsize=8.5, fontweight="bold", color=color) + + fig.subplots_adjust(left=0.08, right=0.86, top=0.82, bottom=0.20, wspace=0.40) + if last_im is not None: + cax = fig.add_axes([0.90, 0.23, 0.018, 0.58]) + cbar = fig.colorbar(last_im, cax=cax) + cbar.set_label("Bootstrap probability") + + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_prefix.with_suffix(ext), bbox_inches="tight") + plt.close(fig) + + +def _balanced_tensor( + matrix_rows: list[dict[str, Any]], + metric: str, +) -> tuple[np.ndarray, list[str], list[str], list[str]] | None: + ok_rows = [r for r in matrix_rows if str(r.get("status", "")).strip() == "ok"] + families = sorted({str(r.get("family", "")).strip() for r in ok_rows if str(r.get("family", "")).strip()}) + decoders = sorted({str(r.get("decoder", "")).strip() for r in ok_rows if str(r.get("decoder", "")).strip()}) + sources = sorted({str(r.get("dataset", "")).strip() for r in ok_rows if str(r.get("dataset", "")).strip()}) + if not families or not decoders or not sources: + return None + + idx = { + (str(r.get("family", "")).strip(), str(r.get("decoder", "")).strip(), str(r.get("dataset", "")).strip()): _f(r.get(metric)) + for r in ok_rows + } + tensor = np.full((len(families), len(decoders), len(sources)), np.nan, dtype=float) + for f_idx, family in enumerate(families): + for d_idx, decoder in enumerate(decoders): + for s_idx, source in enumerate(sources): + tensor[f_idx, d_idx, s_idx] = idx.get((family, decoder, source), float("nan")) + if not np.isfinite(tensor).all(): + return None + return tensor, families, decoders, sources + + +def _variance_decomposition_rows(matrix_rows: list[dict[str, Any]]) -> list[dict[str, Any]]: + metric_labels = [ + ("avg_flip_count", "Average flip count"), + ("nonempty_flip_rate", "Nonempty flip rate"), + ] + component_order = [ + ("family", "Family"), + ("decoder", "Decoder"), + ("source_stack", "Source stack"), + ("family_x_decoder", "Family x decoder"), + ("family_x_source", "Family x source"), + ("decoder_x_source", "Decoder x source"), + ("residual_interaction", "Residual interaction"), + ] + + rows: list[dict[str, Any]] = [] + for metric, metric_label in metric_labels: + packed = _balanced_tensor(matrix_rows, metric) + if packed is None: + continue + tensor, families, decoders, sources = packed + n_family, n_decoder, n_source = tensor.shape + grand = float(np.mean(tensor)) + family_mean = np.mean(tensor, axis=(1, 2)) + decoder_mean = np.mean(tensor, axis=(0, 2)) + source_mean = np.mean(tensor, axis=(0, 1)) + family_decoder_mean = np.mean(tensor, axis=2) + family_source_mean = np.mean(tensor, axis=1) + decoder_source_mean = np.mean(tensor, axis=0) + + sums = { + "family": float(n_decoder * n_source * np.sum((family_mean - grand) ** 2)), + "decoder": float(n_family * n_source * np.sum((decoder_mean - grand) ** 2)), + "source_stack": float(n_family * n_decoder * np.sum((source_mean - grand) ** 2)), + "family_x_decoder": float( + n_source * np.sum((family_decoder_mean - family_mean[:, None] - decoder_mean[None, :] + grand) ** 2) + ), + "family_x_source": float( + n_decoder * np.sum((family_source_mean - family_mean[:, None] - source_mean[None, :] + grand) ** 2) + ), + "decoder_x_source": float( + n_family * np.sum((decoder_source_mean - decoder_mean[:, None] - source_mean[None, :] + grand) ** 2) + ), + } + total = float(np.sum((tensor - grand) ** 2)) + sums["residual_interaction"] = max(0.0, total - sum(sums.values())) + dfs = { + "family": max(0, n_family - 1), + "decoder": max(0, n_decoder - 1), + "source_stack": max(0, n_source - 1), + "family_x_decoder": max(0, (n_family - 1) * (n_decoder - 1)), + "family_x_source": max(0, (n_family - 1) * (n_source - 1)), + "decoder_x_source": max(0, (n_decoder - 1) * (n_source - 1)), + "residual_interaction": max(0, (n_family - 1) * (n_decoder - 1) * (n_source - 1)), + } + + for component_key, component_label in component_order: + ss = sums[component_key] + df = dfs[component_key] + rows.append( + { + "metric": metric, + "metric_label": metric_label, + "component": component_key, + "component_label": component_label, + "degrees_of_freedom": df, + "sum_squares": ss, + "mean_square": ss / df if df else float("nan"), + "variance_share": ss / total if total > 0.0 else float("nan"), + } + ) + + return rows + + +def _plot_variance_decomposition(rows: list[dict[str, Any]], out_prefix: Path) -> None: + try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt # type: ignore + from matplotlib.patches import Patch # type: ignore + except Exception as exc: + print(f"Warning: matplotlib unavailable; skipping variance decomposition plot ({exc}).") + return + + if not rows: + return + + metrics: list[str] = [] + for row in rows: + metric_label = str(row.get("metric_label", "")) + if metric_label and metric_label not in metrics: + metrics.append(metric_label) + component_labels: list[str] = [] + for row in rows: + label = str(row.get("component_label", "")) + if label and label not in component_labels: + component_labels.append(label) + + colors = { + "Family": "#2563EB", + "Decoder": "#DC2626", + "Source stack": "#059669", + "Family x decoder": "#7C3AED", + "Family x source": "#EA580C", + "Decoder x source": "#0891B2", + "Residual interaction": "#64748B", + } + fig, ax = plt.subplots(figsize=(9.4, 3.7), dpi=320, constrained_layout=False) + y_positions = np.arange(len(metrics), dtype=float) + + for y_idx, metric_label in enumerate(metrics): + left = 0.0 + metric_rows = [r for r in rows if str(r.get("metric_label", "")) == metric_label] + for component_label in component_labels: + row = next((r for r in metric_rows if str(r.get("component_label", "")) == component_label), None) + share = _f((row or {}).get("variance_share")) + if not np.isfinite(share): + share = 0.0 + ax.barh( + y_positions[y_idx], + share, + left=left, + height=0.48, + color=colors.get(component_label, "#94A3B8"), + edgecolor="white", + linewidth=0.8, + ) + if share >= 0.075: + ax.text( + left + share / 2.0, + y_positions[y_idx], + f"{share * 100:.0f}%", + ha="center", + va="center", + fontsize=8, + fontweight="bold", + color="white", + ) + left += share + + ax.set_yticks(y_positions) + ax.set_yticklabels(metrics) + ax.set_xlim(0.0, 1.0) + ax.set_xlabel("Share of balanced matrix sum of squares") + ax.set_title("Variance Attribution Across Family, Decoder, and Source Stack", fontsize=11.5, pad=10) + ax.grid(axis="x", alpha=0.22) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.invert_yaxis() + + handles = [Patch(facecolor=colors.get(label, "#94A3B8"), label=label) for label in component_labels] + fig.legend(handles=handles, loc="lower center", ncol=4, frameon=False, bbox_to_anchor=(0.5, -0.04), fontsize=8) + fig.subplots_adjust(left=0.18, right=0.99, top=0.82, bottom=0.30) + + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_prefix.with_suffix(ext), bbox_inches="tight") + plt.close(fig) + + +def _plot_logical_error_rate( + ler_rows: list[dict[str, Any]], + families: list[str], + decoders: list[str], + out_prefix: Path, +) -> None: + try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt # type: ignore + except Exception as exc: + print(f"Warning: matplotlib unavailable; skipping LER plot ({exc}).") + return + + ok_rows = [r for r in ler_rows if str(r.get("status", "")).strip() == "ok"] + if not ok_rows: + return + + cols = max(1, len(families)) + fig, axes = plt.subplots(1, cols, figsize=(5.7 * cols, 4.2), dpi=320, constrained_layout=False, sharey=True) + axes_arr = np.atleast_1d(axes) + palette = {"bp": "#1f77b4", "mwpm": "#2ca02c", "uf": "#d62728"} + + for i, fam in enumerate(families): + ax = axes_arr[i] + fam_rows = [r for r in ok_rows if str(r.get("family", "")) == fam] + datasets = sorted({str(r.get("dataset", "")) for r in fam_rows if str(r.get("dataset", ""))}) + x = np.arange(len(datasets), dtype=float) + width = 0.22 if len(decoders) > 1 else 0.45 + offsets = np.linspace(-width * (len(decoders) - 1), width * (len(decoders) - 1), len(decoders)) if decoders else [] + + for d_idx, decoder in enumerate(decoders): + vals: list[float] = [] + lows: list[float] = [] + highs: list[float] = [] + for dataset in datasets: + row = next( + ( + r + for r in fam_rows + if str(r.get("dataset", "")) == dataset and str(r.get("decoder", "")) == decoder + ), + None, + ) + v = _f((row or {}).get("logical_error_rate")) + lo = _f((row or {}).get("logical_error_ci95_low")) + hi = _f((row or {}).get("logical_error_ci95_high")) + vals.append(v) + lows.append(max(0.0, v - lo) if np.isfinite(v) and np.isfinite(lo) else 0.0) + highs.append(max(0.0, hi - v) if np.isfinite(v) and np.isfinite(hi) else 0.0) + pos = x + float(offsets[d_idx] if len(decoders) > 1 else 0.0) + ax.bar(pos, vals, width=width, color=palette.get(decoder, "#6B7280"), alpha=0.86, label=decoder.upper()) + ax.errorbar(pos, vals, yerr=[lows, highs], fmt="none", ecolor="#111827", elinewidth=0.8, capsize=2.2) + + ax.set_xticks(x) + ax.set_xticklabels(datasets, rotation=22, ha="right", fontsize=8) + ax.set_ylim(0.0, 1.0) + ax.set_title(f"{fam.upper()} logical-parity errors", fontsize=11.5, pad=10) + ax.set_xlabel("Source stack") + if i == 0: + ax.set_ylabel("Logical error rate") + ax.grid(axis="y", alpha=0.25) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.text(-0.12, 1.07, chr(ord("a") + i), transform=ax.transAxes, fontsize=13, fontweight="bold", va="top") + + handles, labels = axes_arr[0].get_legend_handles_labels() + if handles: + fig.legend(handles, labels, loc="lower center", ncol=min(3, len(labels)), frameon=False, bbox_to_anchor=(0.5, -0.02), fontsize=8) + fig.subplots_adjust(left=0.08, right=0.99, top=0.84, bottom=0.30, wspace=0.22) + + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_prefix.with_suffix(ext), bbox_inches="tight") + plt.close(fig) + + def main() -> int: args = parse_args() out_dir = Path(args.out_dir) @@ -217,6 +799,7 @@ def main() -> int: delta_rows: list[dict[str, Any]] = [] combined_matrix_rows: list[dict[str, Any]] = [] combined_source_rows: list[dict[str, Any]] = [] + combined_ler_rows: list[dict[str, Any]] = [] for row in manifest: family = str(row.get("family", "")).strip() @@ -224,11 +807,16 @@ def main() -> int: continue matrix_csv = Path(str(row.get("matrix_csv", "")).strip()) delta_csv = Path(str(row.get("delta_csv", "")).strip()) + results_base_raw = str(row.get("results_base", "")).strip() + ler_csv = Path(results_base_raw) / "03_analysis" / "table_logical_error_rate.csv" if results_base_raw else None if not matrix_csv.exists(): continue matrix_rows = _read_csv(matrix_csv) for mrow in matrix_rows: combined_matrix_rows.append({"family": family, **mrow}) + if ler_csv is not None and ler_csv.exists(): + for ler_row in _read_csv(ler_csv): + combined_ler_rows.append({"family": family, **ler_row}) ok_rows = [r for r in matrix_rows if (r.get("status") or "").strip() == "ok"] decoders = sorted({(r.get("decoder") or "").strip() for r in ok_rows if (r.get("decoder") or "").strip()}) datasets = sorted({(r.get("dataset") or "").strip() for r in ok_rows if (r.get("dataset") or "").strip()}) @@ -333,8 +921,28 @@ def main() -> int: key = (str(row["family"]), str(row["decoder"])) vals = delta_abs_map.get(key, []) row["mean_abs_source_delta_flip"] = float(np.mean(vals)) if vals else float("nan") + ler_src = [ + _f(r.get("logical_error_rate")) + for r in combined_ler_rows + if str(r.get("family", "")) == str(row["family"]) + and str(r.get("decoder", "")) == str(row["decoder"]) + and str(r.get("dataset", "")) != "lidmas_reference" + and str(r.get("status", "")) == "ok" + ] + ler_ref = [ + _f(r.get("logical_error_rate")) + for r in combined_ler_rows + if str(r.get("family", "")) == str(row["family"]) + and str(r.get("decoder", "")) == str(row["decoder"]) + and str(r.get("dataset", "")) == "lidmas_reference" + and str(r.get("status", "")) == "ok" + ] + row["mean_logical_error_rate_sources"] = float(np.mean(ler_src)) if ler_src else float("nan") + row["mean_logical_error_rate_reference"] = float(np.mean(ler_ref)) if ler_ref else float("nan") - families = sorted({str(r["family"]) for r in summary_rows}) + observed_families = {str(r["family"]) for r in summary_rows} + families = [fam for fam in family_order if fam in observed_families] + families.extend(sorted(observed_families - set(families))) for fam in families: fam_rows = [r for r in summary_rows if str(r["family"]) == fam] fam_rows = sorted(fam_rows, key=lambda r: (_f(r["mean_avg_flip_sources"]), str(r["decoder"]))) @@ -382,6 +990,15 @@ def main() -> int: } ) + rank_rows, rank_heatmaps, rank_columns = _rank_stability_by_family( + combined_matrix_rows, + families, + decoders_all, + bootstrap=int(args.rank_bootstrap), + seed=int(args.rank_seed), + ) + variance_rows = _variance_decomposition_rows(combined_matrix_rows) + _write_csv( combined_matrix_rows, [ @@ -441,6 +1058,8 @@ def main() -> int: "mean_avg_flip_reference", "mean_warning_rate_sources", "mean_nonempty_flip_rate_sources", + "mean_logical_error_rate_sources", + "mean_logical_error_rate_reference", "mean_abs_source_delta_flip", "flip_rank_within_family", "norm_flip", @@ -462,16 +1081,59 @@ def main() -> int: ], out_dir / "table_family_delta_effects.csv", ) + _write_csv( + combined_ler_rows, + [ + "family", + "dataset", + "decoder", + "status", + "truth_lines", + "response_lines", + "valid_lines", + "logical_error_count", + "logical_error_rate", + "logical_error_ci95_low", + "logical_error_ci95_high", + "logical_observable", + "truth_model", + "truth_file", + "response_file", + ], + out_dir / "table_logical_error_rate.csv", + ) _write_csv( norm_rows, ["family", "decoder", "norm_flip", "norm_warning", "norm_stack_delta"], out_dir / "table_cross_family_normalized.csv", ) + _write_csv( + rank_rows, + ["family", "decoder", "rank", "rank_prob", "bootstrap_samples", "source_count"], + out_dir / "table_rank_stability_family.csv", + ) + _write_csv( + variance_rows, + [ + "metric", + "metric_label", + "component", + "component_label", + "degrees_of_freedom", + "sum_squares", + "mean_square", + "variance_share", + ], + out_dir / "table_variance_decomposition.csv", + ) _plot_family_tradeoff(summary_rows, families, out_dir / "figure_family_tradeoff") - _plot_family_tradeoff(summary_rows, families, out_dir / "figure_source_vs_lidmas") + _plot_source_vs_lidmas(combined_source_rows, families, out_dir / "figure_source_vs_lidmas") _plot_family_delta_forest(delta_rows, families, out_dir / "figure_family_delta_forest") _plot_normalized_trends(norm_rows, family_order if family_order else families, decoders_all, out_dir / "figure_cross_family_normalized_trends") + _plot_rank_stability_by_family(rank_heatmaps, families, decoders_all, rank_columns, out_dir / "figure_rank_stability_family") + _plot_variance_decomposition(variance_rows, out_dir / "figure_variance_decomposition") + _plot_logical_error_rate(combined_ler_rows, families, decoders_all, out_dir / "figure_logical_error_rate_family") summary_md = out_dir / "summary_code_family_comparison.md" with summary_md.open("w", encoding="utf-8") as f: @@ -479,7 +1141,10 @@ def main() -> int: f.write("This analysis reports:\n\n") f.write("1. within-family decoder tradeoffs (surface and gkp separately),\n") f.write("2. source-vs-reference effect sizes within each family,\n") - f.write("3. cross-family normalized trends (no raw threshold equivalence claims).\n") + f.write("3. cross-family normalized trends (no raw threshold equivalence claims),\n") + f.write("4. source-bootstrap rank stability within each family,\n") + f.write("5. balanced variance attribution across family, decoder, and source stack,\n") + f.write("6. outer-code logical-parity error rates from hidden truth sidecars.\n") return 0 diff --git a/examples/paper_runs/paper_04/scripts/analyze_logical_error_rate.py b/examples/paper_runs/paper_04/scripts/analyze_logical_error_rate.py new file mode 100644 index 0000000..16e1469 --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/analyze_logical_error_rate.py @@ -0,0 +1,296 @@ +#!/usr/bin/env python3 +"""Compute outer-code logical-parity error rates for paper_04 replay outputs.""" + +from __future__ import annotations + +import argparse +import csv +import json +import math +from pathlib import Path +from typing import Any + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--truth-dir", required=True, help="Directory containing truth_*.ndjson sidecars.") + parser.add_argument("--replay-manifest", required=True, help="Replay manifest CSV.") + parser.add_argument("--responses-dir", required=True, help="Directory containing decoder responses.") + parser.add_argument("--out-csv", required=True, help="Output CSV path.") + parser.add_argument("--out-md", required=True, help="Output Markdown path.") + parser.add_argument("--out-prefix", required=True, help="Output figure prefix.") + return parser.parse_args() + + +def _read_csv(path: Path) -> list[dict[str, str]]: + with path.open("r", encoding="utf-8", newline="") as f: + return list(csv.DictReader(f)) + + +def _load_ndjson(path: Path) -> tuple[list[dict[str, Any]], int]: + rows: list[dict[str, Any]] = [] + parse_errors = 0 + with path.open("r", encoding="utf-8") as f: + for line in f: + if not line.strip(): + continue + try: + obj = json.loads(line) + except json.JSONDecodeError: + parse_errors += 1 + continue + rows.append(obj if isinstance(obj, dict) else {}) + return rows, parse_errors + + +def _safe_int(value: Any) -> int: + try: + return int(value) + except Exception: + return 0 + + +def _correction_indices(response: dict[str, Any]) -> set[int]: + correction = response.get("correction", {}) + if not isinstance(correction, dict): + correction = {} + flips = correction.get("qubit_flips_x") + if not isinstance(flips, list): + flips = correction.get("qubit_flips", []) + if not isinstance(flips, list): + flips = [] + return {_safe_int(item) for item in flips} + + +def _parity(indices: set[int], logical_indices: list[int]) -> int: + value = 0 + for idx in logical_indices: + if idx in indices: + value ^= 1 + return value + + +def _wilson(k: int, n: int, z: float = 1.959963984540054) -> tuple[float, float]: + if n <= 0: + return float("nan"), float("nan") + phat = k / n + denom = 1.0 + z * z / n + center = (phat + z * z / (2.0 * n)) / denom + half = z * math.sqrt((phat * (1.0 - phat) / n) + (z * z / (4.0 * n * n))) / denom + return max(0.0, center - half), min(1.0, center + half) + + +def _status(*, truth_path: Path, response_path: Path, truth_parse_errors: int, response_parse_errors: int, truth_lines: int, response_lines: int) -> str: + if not truth_path.exists(): + return "missing_truth" + if not response_path.exists(): + return "missing_response" + if truth_parse_errors or response_parse_errors: + return "parse_errors" + if truth_lines != response_lines: + return "line_count_mismatch" + return "ok" + + +def _analyze_cell(dataset: str, decoder: str, truth_path: Path, response_path: Path) -> dict[str, Any]: + if not truth_path.exists() or not response_path.exists(): + return { + "dataset": dataset, + "decoder": decoder, + "status": _status( + truth_path=truth_path, + response_path=response_path, + truth_parse_errors=0, + response_parse_errors=0, + truth_lines=0, + response_lines=0, + ), + "truth_lines": 0, + "response_lines": 0, + "valid_lines": 0, + "logical_error_count": 0, + "logical_error_rate": float("nan"), + "logical_error_ci95_low": float("nan"), + "logical_error_ci95_high": float("nan"), + "logical_observable": "", + "truth_model": "", + "truth_file": truth_path.name, + "response_file": response_path.name, + } + + truth_rows, truth_parse_errors = _load_ndjson(truth_path) + response_rows, response_parse_errors = _load_ndjson(response_path) + paired = min(len(truth_rows), len(response_rows)) + logical_errors = 0 + valid = 0 + logical_observable = "" + truth_model = "" + + for idx in range(paired): + truth = truth_rows[idx] + response = response_rows[idx] + logical_indices = truth.get("logical_indices", []) + if not isinstance(logical_indices, list): + continue + logical_indices_int = [_safe_int(item) for item in logical_indices] + truth_value = _safe_int(truth.get("logical_truth", 0)) & 1 + correction_value = _parity(_correction_indices(response), logical_indices_int) + residual_value = truth_value ^ correction_value + logical_errors += residual_value + valid += 1 + logical_observable = str(truth.get("logical_observable", logical_observable)) + truth_model = str(truth.get("truth_model", truth_model)) + + low, high = _wilson(logical_errors, valid) + return { + "dataset": dataset, + "decoder": decoder, + "status": _status( + truth_path=truth_path, + response_path=response_path, + truth_parse_errors=truth_parse_errors, + response_parse_errors=response_parse_errors, + truth_lines=len(truth_rows), + response_lines=len(response_rows), + ), + "truth_lines": len(truth_rows), + "response_lines": len(response_rows), + "valid_lines": valid, + "logical_error_count": logical_errors, + "logical_error_rate": logical_errors / valid if valid else float("nan"), + "logical_error_ci95_low": low, + "logical_error_ci95_high": high, + "logical_observable": logical_observable, + "truth_model": truth_model, + "truth_file": truth_path.name, + "response_file": response_path.name, + } + + +def _fmt(value: Any) -> Any: + if isinstance(value, float): + if not math.isfinite(value): + return "" + return f"{value:.6f}" + return value + + +def _write_csv(rows: list[dict[str, Any]], out_csv: Path) -> None: + fields = [ + "dataset", + "decoder", + "status", + "truth_lines", + "response_lines", + "valid_lines", + "logical_error_count", + "logical_error_rate", + "logical_error_ci95_low", + "logical_error_ci95_high", + "logical_observable", + "truth_model", + "truth_file", + "response_file", + ] + out_csv.parent.mkdir(parents=True, exist_ok=True) + with out_csv.open("w", encoding="utf-8", newline="") as f: + writer = csv.DictWriter(f, fieldnames=fields) + writer.writeheader() + for row in rows: + writer.writerow({field: _fmt(row.get(field, "")) for field in fields}) + + +def _write_md(rows: list[dict[str, Any]], out_md: Path) -> None: + fields = ["dataset", "decoder", "status", "valid_lines", "logical_error_count", "logical_error_rate", "logical_error_ci95_low", "logical_error_ci95_high"] + with out_md.open("w", encoding="utf-8") as f: + f.write("| " + " | ".join(fields) + " |\n") + f.write("| " + " | ".join(["---"] * len(fields)) + " |\n") + for row in rows: + f.write("| " + " | ".join(str(_fmt(row.get(field, ""))) for field in fields) + " |\n") + + +def _plot(rows: list[dict[str, Any]], out_prefix: Path) -> None: + try: + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt # type: ignore + import numpy as np # type: ignore + except Exception as exc: + print(f"Warning: matplotlib unavailable; skipping LER plot ({exc}).") + return + + ok_rows = [r for r in rows if str(r.get("status", "")) == "ok"] + if not ok_rows: + return + + datasets = sorted({str(r.get("dataset", "")) for r in ok_rows}) + decoders = sorted({str(r.get("decoder", "")) for r in ok_rows}) + colors = {"bp": "#1f77b4", "mwpm": "#2ca02c", "uf": "#d62728"} + + x = np.arange(len(datasets), dtype=float) + width = 0.22 if len(decoders) > 1 else 0.45 + offsets = np.linspace(-width * (len(decoders) - 1), width * (len(decoders) - 1), len(decoders)) if decoders else [] + fig, ax = plt.subplots(figsize=(7.4, 4.1), dpi=320, constrained_layout=True) + + for d_idx, decoder in enumerate(decoders): + vals = [] + lows = [] + highs = [] + for dataset in datasets: + row = next((r for r in ok_rows if str(r.get("dataset", "")) == dataset and str(r.get("decoder", "")) == decoder), None) + v = float(row.get("logical_error_rate", float("nan"))) if row else float("nan") + lo = float(row.get("logical_error_ci95_low", float("nan"))) if row else float("nan") + hi = float(row.get("logical_error_ci95_high", float("nan"))) if row else float("nan") + vals.append(v) + lows.append(max(0.0, v - lo) if np.isfinite(v) and np.isfinite(lo) else 0.0) + highs.append(max(0.0, hi - v) if np.isfinite(v) and np.isfinite(hi) else 0.0) + pos = x + float(offsets[d_idx] if len(decoders) > 1 else 0.0) + ax.bar(pos, vals, width=width, color=colors.get(decoder, "#6B7280"), label=decoder.upper(), alpha=0.86) + ax.errorbar(pos, vals, yerr=[lows, highs], fmt="none", ecolor="#111827", elinewidth=0.8, capsize=2.4) + + ax.set_xticks(x) + ax.set_xticklabels(datasets, rotation=20, ha="right") + ax.set_ylim(0.0, 1.0) + ax.set_ylabel("Logical error rate") + ax.set_title("Outer-code logical-parity error rate") + ax.grid(axis="y", alpha=0.25) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.legend(frameon=False, ncol=min(3, len(decoders))) + + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_prefix.with_suffix(ext), bbox_inches="tight") + plt.close(fig) + + +def main() -> int: + args = parse_args() + truth_dir = Path(args.truth_dir) + responses_dir = Path(args.responses_dir) + manifest_rows = _read_csv(Path(args.replay_manifest)) + + rows: list[dict[str, Any]] = [] + for row in manifest_rows: + dataset = str(row.get("dataset", "")).strip() + decoder = str(row.get("decoder", "")).strip() + response_file = str(row.get("response_file", "")).strip() + if not dataset or not decoder or not response_file: + continue + rows.append( + _analyze_cell( + dataset=dataset, + decoder=decoder, + truth_path=truth_dir / f"truth_{dataset}.ndjson", + response_path=responses_dir / response_file, + ) + ) + + _write_csv(rows, Path(args.out_csv)) + _write_md(rows, Path(args.out_md)) + _plot(rows, Path(args.out_prefix)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/paper_runs/paper_04/scripts/compose_journal_results_figure.py b/examples/paper_runs/paper_04/scripts/compose_journal_results_figure.py new file mode 100644 index 0000000..240e162 --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/compose_journal_results_figure.py @@ -0,0 +1,565 @@ +#!/usr/bin/env python3 +"""Compose a journal-style multi-panel summary figure for paper_04 results.""" + +from __future__ import annotations + +import argparse +import csv +import math +import os +import shutil +import tempfile +from pathlib import Path +from typing import Any + +import numpy as np + + +DECODER_ORDER = ["bp", "mwpm", "uf"] +FAMILY_ORDER = ["surface", "gkp"] +SOURCE_ORDER = ["cirq", "pennylane", "qiskit"] +DECODER_COLORS = { + "bp": "#2563EB", + "mwpm": "#059669", + "uf": "#DC2626", +} +SOURCE_COLORS = { + "cirq": "#7C3AED", + "pennylane": "#0284C7", + "qiskit": "#16A34A", +} +FAMILY_MARKERS = { + "surface": "o", + "gkp": "s", +} + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--analysis-dir", required=True, help="Directory containing paper_04 analysis CSV files.") + parser.add_argument("--out-prefix", required=True, help="Output prefix without extension.") + parser.add_argument("--manuscript-dir", help="Optional manuscript figure directory to receive copied outputs.") + parser.add_argument("--write-standalone", action="store_true", help="Overwrite standalone manuscript result figures using the journal visual style.") + return parser.parse_args() + + +def read_csv(path: Path) -> list[dict[str, str]]: + with path.open("r", encoding="utf-8", newline="") as f: + return list(csv.DictReader(f)) + + +def f(row: dict[str, Any], key: str) -> float: + try: + return float(row.get(key, "nan")) + except (TypeError, ValueError): + return float("nan") + + +def wilson(k: int, n: int, z: float = 1.959963984540054) -> tuple[float, float, float]: + if n <= 0: + return float("nan"), float("nan"), float("nan") + phat = k / n + denom = 1.0 + z * z / n + center = (phat + z * z / (2.0 * n)) / denom + half = z * math.sqrt((phat * (1.0 - phat) + z * z / (4.0 * n)) / n) / denom + return phat, max(0.0, center - half), min(1.0, center + half) + + +def panel_label(ax: Any, label: str, x: float = -0.13, y: float = 1.08) -> None: + ax.text(x, y, label, transform=ax.transAxes, fontsize=8.6, fontweight="bold", va="top", ha="left") + + +def tidy_axes(ax: Any, *, grid: str | None = "y") -> None: + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + if grid: + ax.grid(axis=grid, color="#E5E7EB", linewidth=0.45, zorder=0) + ax.tick_params(width=0.6, length=2.5, color="#374151") + + +def save_outputs(fig: Any, out_prefix: Path, manuscript_dir: Path | None = None) -> list[Path]: + out_files: list[Path] = [] + out_prefix.parent.mkdir(parents=True, exist_ok=True) + for ext in (".pdf", ".png", ".svg"): + out = out_prefix.with_suffix(ext) + fig.savefig(out, bbox_inches="tight") + out_files.append(out) + if manuscript_dir is not None: + manuscript_dir.mkdir(parents=True, exist_ok=True) + for out in out_files: + shutil.copy2(out, manuscript_dir / out.name) + return out_files + + +def source_mean_ler(ler_rows: list[dict[str, str]]) -> dict[tuple[str, str], tuple[float, float, float]]: + out: dict[tuple[str, str], tuple[float, float, float]] = {} + for fam in FAMILY_ORDER: + for dec in DECODER_ORDER: + rows = [ + r + for r in ler_rows + if r.get("family") == fam and r.get("decoder") == dec and r.get("dataset") != "lidmas_reference" + ] + k = sum(int(float(r.get("logical_error_count", "0"))) for r in rows) + n = sum(int(float(r.get("valid_lines", "0"))) for r in rows) + out[(fam, dec)] = wilson(k, n) + return out + + +def add_panel_a_tradeoff(ax: Any, summary_rows: list[dict[str, str]], label: str | None = "a") -> None: + for fam in FAMILY_ORDER: + rows = [r for r in summary_rows if r.get("family") == fam] + rows = sorted(rows, key=lambda r: DECODER_ORDER.index(r.get("decoder", ""))) + xs = [f(r, "mean_avg_flip_sources") for r in rows] + ys = [f(r, "mean_logical_error_rate_sources") for r in rows] + ax.plot(xs, ys, color="#9CA3AF", linewidth=0.8, zorder=1) + for r, x, y in zip(rows, xs, ys): + dec = r.get("decoder", "") + ax.scatter( + x, + y, + s=42, + marker=FAMILY_MARKERS[fam], + facecolor=DECODER_COLORS.get(dec, "#6B7280"), + edgecolor="white", + linewidth=0.65, + zorder=3, + ) + ax.text(x + 0.05, y + 0.002, dec.upper(), fontsize=5.5, color="#111827", va="bottom") + + ax.set_xlabel("mean flips") + ax.set_ylabel("logical-parity error rate") + ax.set_ylim(0.425, 0.515) + ax.set_xlim(1.25, 6.15) + if label: + panel_label(ax, label) + tidy_axes(ax, grid="both") + + from matplotlib.lines import Line2D # type: ignore + + handles = [ + Line2D([0], [0], marker=FAMILY_MARKERS["surface"], color="none", markerfacecolor="#6B7280", markeredgecolor="white", markersize=5.0, label="surface"), + Line2D([0], [0], marker=FAMILY_MARKERS["gkp"], color="none", markerfacecolor="#6B7280", markeredgecolor="white", markersize=5.0, label="GKP"), + ] + ax.legend(handles=handles, frameon=False, loc="lower right", fontsize=5.2, handlelength=0.8, borderpad=0.2) + + +def add_delta_axis(ax: Any, rows: list[dict[str, str]], fam: str, show_ylabel: bool) -> None: + fam_rows = [r for r in rows if r.get("family") == fam] + fam_rows = sorted( + fam_rows, + key=lambda r: (DECODER_ORDER.index(r.get("decoder", "")), SOURCE_ORDER.index(r.get("source_dataset", ""))), + ) + y = np.arange(len(fam_rows), dtype=float) + for yi, r in zip(y, fam_rows): + src = r.get("source_dataset", "") + dec = r.get("decoder", "") + mean = f(r, "delta_mean_source_minus_reference") + lo = f(r, "delta_ci95_low") + hi = f(r, "delta_ci95_high") + ax.plot([lo, hi], [yi, yi], color=SOURCE_COLORS.get(src, "#6B7280"), linewidth=0.85, alpha=0.95) + ax.scatter(mean, yi, s=13, color=SOURCE_COLORS.get(src, "#6B7280"), edgecolor="white", linewidth=0.35, zorder=3) + ax.axvline(0.0, color="#111827", linewidth=0.6, alpha=0.7) + ax.set_title(fam.upper(), fontsize=6.2, pad=2) + ax.set_yticks(y) + if show_ylabel: + labels = [f"{r.get('decoder','').upper()}-{r.get('source_dataset','')[:1].upper()}" for r in fam_rows] + ax.set_yticklabels(labels, fontsize=4.8) + else: + ax.set_yticklabels([]) + ax.invert_yaxis() + ax.set_xlabel(r"$\Delta$ flips") + tidy_axes(ax, grid="x") + if fam == "surface": + ax.set_xlim(-0.26, 0.31) + else: + ax.set_xlim(-2.45, 0.72) + + +def add_panel_c_ler(ax: Any, ler_rows: list[dict[str, str]], label: str | None = "c") -> None: + agg = source_mean_ler(ler_rows) + family_x = {"surface": 0.0, "gkp": 1.0} + offsets = {"bp": -0.17, "mwpm": 0.0, "uf": 0.17} + for fam in FAMILY_ORDER: + for dec in DECODER_ORDER: + rate, lo, hi = agg[(fam, dec)] + x = family_x[fam] + offsets[dec] + ax.errorbar( + [x], + [rate], + yerr=[[rate - lo], [hi - rate]], + fmt=FAMILY_MARKERS[fam], + markersize=4.8, + color=DECODER_COLORS[dec], + markeredgecolor="white", + markeredgewidth=0.55, + elinewidth=0.8, + capsize=2.2, + zorder=3, + ) + ax.set_xticks([0, 1]) + ax.set_xticklabels(["surface", "GKP"]) + ax.set_ylabel(r"$\lambda$ (source aggregate)") + ax.set_ylim(0.39, 0.53) + if label: + panel_label(ax, label) + tidy_axes(ax, grid="y") + + from matplotlib.lines import Line2D # type: ignore + + handles = [ + Line2D([0], [0], marker="o", color=DECODER_COLORS[d], label=d.upper(), markersize=4.2, linewidth=0.9) + for d in DECODER_ORDER + ] + ax.legend(handles=handles, frameon=False, loc="upper left", fontsize=5.2, ncol=3, handlelength=1.0, columnspacing=0.6) + + +def add_panel_d_rank(ax: Any, rank_rows: list[dict[str, str]], label: str | None = "d") -> Any: + rows: list[tuple[str, str]] = [(fam, dec) for fam in FAMILY_ORDER for dec in DECODER_ORDER] + mat = np.zeros((len(rows), 3), dtype=float) + for i, (fam, dec) in enumerate(rows): + for r in rank_rows: + if r.get("family") == fam and r.get("decoder") == dec: + rank = int(float(r.get("rank", "0"))) + if 1 <= rank <= 3: + mat[i, rank - 1] = f(r, "rank_prob") + im = ax.imshow(mat, cmap="Blues", vmin=0.0, vmax=1.0, aspect="auto") + for i in range(mat.shape[0]): + for j in range(mat.shape[1]): + val = mat[i, j] + color = "white" if val > 0.55 else "#111827" + ax.text(j, i, f"{val:.2f}", ha="center", va="center", fontsize=4.9, color=color, fontweight="bold" if val > 0.55 else "normal") + ax.set_xticks([0, 1, 2]) + ax.set_xticklabels(["1", "2", "3"]) + ax.set_yticks(np.arange(len(rows))) + ax.set_yticklabels([f"{fam[0].upper()} {dec.upper()}" for fam, dec in rows], fontsize=5.1) + ax.set_xlabel("rank") + ax.set_ylabel("family decoder") + if label: + panel_label(ax, label) + for spine in ax.spines.values(): + spine.set_visible(False) + ax.tick_params(length=0) + return im + + +def add_panel_e_variance(ax: Any, variance_rows: list[dict[str, str]], label: str | None = "e") -> None: + metrics = ["avg_flip_count", "nonempty_flip_rate"] + metric_labels = {"avg_flip_count": "mean flips", "nonempty_flip_rate": "nonempty rate"} + components = [ + "family", + "decoder", + "source_stack", + "family_x_decoder", + "family_x_source", + "decoder_x_source", + "residual_interaction", + ] + comp_labels = { + "family": "family", + "decoder": "decoder", + "source_stack": "source", + "family_x_decoder": "fam x dec", + "family_x_source": "fam x src", + "decoder_x_source": "dec x src", + "residual_interaction": "resid.", + } + comp_colors = { + "family": "#2563EB", + "decoder": "#DC2626", + "source_stack": "#059669", + "family_x_decoder": "#8B5CF6", + "family_x_source": "#F97316", + "decoder_x_source": "#0891B2", + "residual_interaction": "#64748B", + } + y = np.arange(len(metrics), dtype=float) + for yi, metric in zip(y, metrics): + left = 0.0 + for comp in components: + row = next((r for r in variance_rows if r.get("metric") == metric and r.get("component") == comp), None) + share = f(row or {}, "variance_share") + if not np.isfinite(share): + share = 0.0 + ax.barh(yi, share, left=left, height=0.48, color=comp_colors[comp], edgecolor="white", linewidth=0.35) + if share >= 0.10: + ax.text(left + share / 2, yi, f"{100*share:.0f}%", ha="center", va="center", fontsize=4.9, color="white", fontweight="bold") + left += share + ax.set_xlim(0.0, 1.0) + ax.set_yticks(y) + ax.set_yticklabels([metric_labels[m] for m in metrics]) + ax.set_xlabel("variance share") + ax.invert_yaxis() + if label: + panel_label(ax, label) + tidy_axes(ax, grid="x") + + from matplotlib.patches import Patch # type: ignore + + handles = [Patch(facecolor=comp_colors[c], label=comp_labels[c]) for c in components] + ax.legend(handles=handles, frameon=False, loc="lower center", bbox_to_anchor=(0.5, -0.54), ncol=4, fontsize=4.8, handlelength=1.0, columnspacing=0.8) + + +def add_panel_f_normalized(ax: Any, norm_rows: list[dict[str, str]], label: str | None = "f") -> None: + metrics = [("norm_flip", "flip"), ("norm_stack_delta", "source delta")] + rows: list[tuple[str, str, str]] = [] + for metric, label in metrics: + for dec in DECODER_ORDER: + rows.append((metric, label, dec)) + y = np.arange(len(rows), dtype=float) + for yi, (metric, _label, dec) in zip(y, rows): + values = {} + for fam in FAMILY_ORDER: + row = next((r for r in norm_rows if r.get("family") == fam and r.get("decoder") == dec), None) + values[fam] = f(row or {}, metric) + ax.plot([values["surface"], values["gkp"]], [yi, yi], color=DECODER_COLORS[dec], linewidth=0.85, alpha=0.55) + ax.scatter(values["surface"], yi, s=17, marker="o", color=DECODER_COLORS[dec], edgecolor="white", linewidth=0.4, zorder=3) + ax.scatter(values["gkp"], yi, s=20, marker="s", color=DECODER_COLORS[dec], edgecolor="white", linewidth=0.4, zorder=3) + ax.set_xlim(-0.05, 1.05) + ax.set_yticks(y) + ax.set_yticklabels([f"{label} {dec.upper()}" for _, label, dec in rows], fontsize=5.1) + ax.set_xlabel("within-family normalized value") + ax.invert_yaxis() + if label: + panel_label(ax, label) + tidy_axes(ax, grid="x") + + from matplotlib.lines import Line2D # type: ignore + + handles = [ + Line2D([0], [0], marker="o", color="none", markerfacecolor="#6B7280", markeredgecolor="white", markersize=4.2, label="surface"), + Line2D([0], [0], marker="s", color="none", markerfacecolor="#6B7280", markeredgecolor="white", markersize=4.2, label="GKP"), + ] + ax.legend(handles=handles, frameon=False, loc="upper right", fontsize=5.2, handlelength=0.9) + + +def add_source_vs_reference_axis(ax: Any, rows: list[dict[str, str]], fam: str, show_ylabel: bool) -> None: + fam_rows = [r for r in rows if r.get("family") == fam] + x_base = {dec: i for i, dec in enumerate(DECODER_ORDER)} + offsets = {"cirq": -0.16, "pennylane": 0.0, "qiskit": 0.16} + ref_seen: set[str] = set() + for dec in DECODER_ORDER: + dec_rows = [r for r in fam_rows if r.get("decoder") == dec] + if not dec_rows: + continue + x = x_base[dec] + ref = f(dec_rows[0], "avg_flip_count_reference") + ax.scatter( + x, + ref, + marker="D", + s=24, + color="#111827", + edgecolor="white", + linewidth=0.45, + zorder=4, + label="LiDMaS+ reference" if not ref_seen else None, + ) + ref_seen.add("reference") + for r in dec_rows: + src = r.get("source_dataset", "") + x_src = x + offsets.get(src, 0.0) + y_src = f(r, "avg_flip_count_source") + ax.plot([x, x_src], [ref, y_src], color=SOURCE_COLORS.get(src, "#6B7280"), linewidth=0.75, alpha=0.55, zorder=1) + ax.scatter( + x_src, + y_src, + marker="o", + s=25, + color=SOURCE_COLORS.get(src, "#6B7280"), + edgecolor="white", + linewidth=0.45, + zorder=3, + label=src if dec == DECODER_ORDER[0] else None, + ) + ax.set_title(fam.upper(), fontsize=8.0, pad=5) + ax.set_xticks([x_base[d] for d in DECODER_ORDER]) + ax.set_xticklabels([d.upper() for d in DECODER_ORDER]) + ax.set_xlabel("decoder") + if show_ylabel: + ax.set_ylabel("mean flips") + else: + ax.set_ylabel("") + if fam == "gkp": + ax.set_ylim(0.55, 4.65) + else: + ax.set_ylim(2.2, 6.25) + tidy_axes(ax, grid="y") + + +def render_standalone_figures( + analysis_dir: Path, + manuscript_dir: Path | None, + summary_rows: list[dict[str, str]], + delta_rows: list[dict[str, str]], + source_rows: list[dict[str, str]], + ler_rows: list[dict[str, str]], + rank_rows: list[dict[str, str]], + variance_rows: list[dict[str, str]], + norm_rows: list[dict[str, str]], +) -> list[Path]: + import matplotlib.pyplot as plt # type: ignore + from matplotlib.lines import Line2D # type: ignore + + written: list[Path] = [] + + fig, axes = plt.subplots(1, 2, figsize=(7.05, 2.75), constrained_layout=False) + add_source_vs_reference_axis(axes[0], source_rows, "gkp", True) + add_source_vs_reference_axis(axes[1], source_rows, "surface", False) + handles = [ + Line2D([0], [0], marker="D", color="none", markerfacecolor="#111827", markeredgecolor="white", markersize=4.3, label="LiDMaS+ reference"), + *[ + Line2D([0], [0], marker="o", color=SOURCE_COLORS[src], markerfacecolor=SOURCE_COLORS[src], markeredgecolor="white", markersize=4.3, linewidth=0.9, label=src) + for src in SOURCE_ORDER + ], + ] + fig.legend(handles=handles, frameon=False, ncol=4, loc="lower center", bbox_to_anchor=(0.5, -0.03), fontsize=6.0, handlelength=1.0, columnspacing=0.9) + fig.subplots_adjust(bottom=0.22) + written.extend(save_outputs(fig, analysis_dir / "figure_source_vs_lidmas", manuscript_dir)) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(4.65, 3.25), constrained_layout=True) + add_panel_a_tradeoff(ax, summary_rows, label=None) + written.extend(save_outputs(fig, analysis_dir / "figure_family_tradeoff", manuscript_dir)) + plt.close(fig) + + fig = plt.figure(figsize=(7.05, 3.0), constrained_layout=True) + gs = fig.add_gridspec(1, 2, wspace=0.18) + ax1 = fig.add_subplot(gs[0, 0]) + ax2 = fig.add_subplot(gs[0, 1]) + add_delta_axis(ax1, delta_rows, "gkp", True) + add_delta_axis(ax2, delta_rows, "surface", False) + handles = [ + Line2D([0], [0], marker="o", color=SOURCE_COLORS[src], markerfacecolor=SOURCE_COLORS[src], markeredgecolor="white", markersize=4.3, linewidth=0.9, label=src) + for src in SOURCE_ORDER + ] + ax2.legend(handles=handles, frameon=False, loc="upper right", fontsize=5.8, handlelength=0.9) + written.extend(save_outputs(fig, analysis_dir / "figure_family_delta_forest", manuscript_dir)) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(4.4, 3.0), constrained_layout=True) + add_panel_f_normalized(ax, norm_rows, label=None) + written.extend(save_outputs(fig, analysis_dir / "figure_cross_family_normalized_trends", manuscript_dir)) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(4.2, 3.1), constrained_layout=True) + add_panel_d_rank(ax, rank_rows, label=None) + written.extend(save_outputs(fig, analysis_dir / "figure_rank_stability_family", manuscript_dir)) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(5.4, 2.75), constrained_layout=True) + add_panel_e_variance(ax, variance_rows, label=None) + written.extend(save_outputs(fig, analysis_dir / "figure_variance_decomposition", manuscript_dir)) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(4.25, 3.0), constrained_layout=True) + add_panel_c_ler(ax, ler_rows, label=None) + written.extend(save_outputs(fig, analysis_dir / "figure_logical_error_rate_family", manuscript_dir)) + plt.close(fig) + + return written + + +def main() -> int: + args = parse_args() + analysis_dir = Path(args.analysis_dir) + out_prefix = Path(args.out_prefix) + out_prefix.parent.mkdir(parents=True, exist_ok=True) + mpl_config_dir = Path(os.environ.get("MPLCONFIGDIR", Path(tempfile.gettempdir()) / "lidmas_matplotlib")) + mpl_config_dir.mkdir(parents=True, exist_ok=True) + os.environ.setdefault("MPLCONFIGDIR", str(mpl_config_dir)) + + summary_rows = read_csv(analysis_dir / "table_family_decoder_summary.csv") + delta_rows = read_csv(analysis_dir / "table_family_delta_effects.csv") + source_rows = read_csv(analysis_dir / "table_source_vs_lidmas.csv") + ler_rows = read_csv(analysis_dir / "table_logical_error_rate.csv") + rank_rows = read_csv(analysis_dir / "table_rank_stability_family.csv") + variance_rows = read_csv(analysis_dir / "table_variance_decomposition.csv") + norm_rows = read_csv(analysis_dir / "table_cross_family_normalized.csv") + + import matplotlib # type: ignore + + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt # type: ignore + + plt.rcParams.update( + { + "pdf.fonttype": 42, + "ps.fonttype": 42, + "font.family": "Arial", + "font.size": 6.2, + "axes.labelsize": 6.2, + "axes.titlesize": 6.2, + "xtick.labelsize": 5.4, + "ytick.labelsize": 5.4, + "legend.fontsize": 5.2, + "axes.linewidth": 0.65, + "savefig.dpi": 420, + } + ) + + fig = plt.figure(figsize=(7.05, 6.35), constrained_layout=False) + gs = fig.add_gridspec(3, 2, height_ratios=[1.05, 1.0, 0.92], hspace=0.62, wspace=0.42) + + ax_a = fig.add_subplot(gs[0, 0]) + add_panel_a_tradeoff(ax_a, summary_rows) + + sub_b = gs[0, 1].subgridspec(1, 2, wspace=0.22) + ax_b1 = fig.add_subplot(sub_b[0, 0]) + ax_b2 = fig.add_subplot(sub_b[0, 1]) + add_delta_axis(ax_b1, delta_rows, "gkp", True) + add_delta_axis(ax_b2, delta_rows, "surface", False) + panel_label(ax_b1, "b", x=-0.22) + from matplotlib.lines import Line2D # type: ignore + + ax_b2.legend( + handles=[ + Line2D([0], [0], marker="o", color=SOURCE_COLORS[s], label=s, markersize=3.2, linewidth=0.8) + for s in SOURCE_ORDER + ], + frameon=False, + loc="upper right", + fontsize=4.8, + handlelength=0.9, + ) + + ax_c = fig.add_subplot(gs[1, 0]) + add_panel_c_ler(ax_c, ler_rows) + + ax_d = fig.add_subplot(gs[1, 1]) + add_panel_d_rank(ax_d, rank_rows) + + ax_e = fig.add_subplot(gs[2, 0]) + add_panel_e_variance(ax_e, variance_rows) + + ax_f = fig.add_subplot(gs[2, 1]) + add_panel_f_normalized(ax_f, norm_rows) + + fig.subplots_adjust(left=0.075, right=0.985, top=0.972, bottom=0.095) + + manuscript_dir = Path(args.manuscript_dir) if args.manuscript_dir else None + out_files = save_outputs(fig, out_prefix, manuscript_dir) + plt.close(fig) + + if args.write_standalone: + out_files.extend( + render_standalone_figures( + analysis_dir, + manuscript_dir, + summary_rows, + delta_rows, + source_rows, + ler_rows, + rank_rows, + variance_rows, + norm_rows, + ) + ) + + for out in out_files: + print(f"Wrote {out}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/paper_runs/paper_04/scripts/compose_study_architecture.py b/examples/paper_runs/paper_04/scripts/compose_study_architecture.py new file mode 100644 index 0000000..4f67591 --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/compose_study_architecture.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python3 +"""Compose two Mermaid-rendered architecture panels into one paper figure.""" + +from __future__ import annotations + +from pathlib import Path + +from PIL import Image, ImageDraw, ImageFont + + +ROOT = Path("examples/paper_runs/paper_04/results/03_analysis") + + +def _font(size: int, *, bold: bool = False) -> ImageFont.FreeTypeFont | ImageFont.ImageFont: + candidates = [ + Path("/System/Library/Fonts/Supplemental/Arial Bold.ttf" if bold else "/System/Library/Fonts/Supplemental/Arial.ttf"), + Path("/System/Library/Fonts/ArialHB.ttc"), + ] + for path in candidates: + if path.exists(): + try: + return ImageFont.truetype(str(path), size=size) + except OSError: + continue + return ImageFont.load_default() + + +def _fit_width(image: Image.Image, target_width: int) -> Image.Image: + if image.width == target_width: + return image + scale = target_width / image.width + return image.resize((target_width, round(image.height * scale)), Image.Resampling.LANCZOS) + + +def _orange_box(image: Image.Image) -> tuple[int, int, int, int]: + """Return the largest orange outline component in a rendered Mermaid panel.""" + pixels = image.load() + orange: set[tuple[int, int]] = set() + for y in range(image.height): + for x in range(image.width): + red, green, blue = pixels[x, y] + if red > 180 and 45 < green < 150 and blue < 90: + orange.add((x, y)) + + best: tuple[int, int, int, int, int] | None = None + while orange: + start = orange.pop() + stack = [start] + xs: list[int] = [] + ys: list[int] = [] + while stack: + x, y = stack.pop() + xs.append(x) + ys.append(y) + for nx, ny in ((x + 1, y), (x - 1, y), (x, y + 1), (x, y - 1)): + if (nx, ny) in orange: + orange.remove((nx, ny)) + stack.append((nx, ny)) + if len(xs) > 20: + component = (len(xs), min(xs), min(ys), max(xs), max(ys)) + if best is None or component[0] > best[0]: + best = component + + if best is None: + raise RuntimeError("No orange contract box found in Mermaid panel.") + _, left, top, right, bottom = best + return left, top, right, bottom + + +def main() -> int: + generation = Image.open(ROOT / "figure_study_architecture_generation.png").convert("RGB") + replay = Image.open(ROOT / "figure_study_architecture_replay.png").convert("RGB") + + margin = 48 + arrow_gap = 96 + target_width = max(generation.width, replay.width, 1000) + generation = _fit_width(generation, target_width) + replay = _fit_width(replay, target_width) + + bridge_font = _font(34, bold=True) + + width = target_width + margin * 2 + height = margin + generation.height + arrow_gap + replay.height + margin + canvas = Image.new("RGB", (width, height), "#FFFFFF") + draw = ImageDraw.Draw(canvas) + + y = margin + generation_y = y + canvas.paste(generation, (margin, generation_y)) + y += generation.height + + label = "contract-preserving replay" + bbox = draw.textbbox((0, 0), label, font=bridge_font) + label_w = bbox[2] - bbox[0] + + replay_y = y + arrow_gap + canvas.paste(replay, (margin, replay_y)) + + top_contract = _orange_box(generation) + bottom_contract = _orange_box(replay) + start_x = margin + (top_contract[0] + top_contract[2]) // 2 + start_y = generation_y + top_contract[3] + 4 + end_x = margin + (bottom_contract[0] + bottom_contract[2]) // 2 + end_y = replay_y + bottom_contract[1] - 2 + connector_y = generation_y + generation.height + 58 + + label_x = end_x + (start_x - end_x - label_w) // 2 + draw.text((label_x, generation_y + generation.height + 8), label, fill="#334155", font=bridge_font) + draw.line((start_x, start_y, start_x, connector_y, end_x, connector_y, end_x, end_y), fill="#334155", width=7) + draw.polygon([(end_x - 20, end_y - 24), (end_x + 20, end_y - 24), (end_x, end_y + 4)], fill="#334155") + y += arrow_gap + + out_png = ROOT / "figure_study_architecture.png" + out_pdf = ROOT / "figure_study_architecture.pdf" + canvas.save(out_png, optimize=True) + canvas.save(out_pdf, "PDF", resolution=300.0) + manuscript = ROOT / "manuscript_figures" + manuscript.mkdir(exist_ok=True) + manuscript_png = manuscript / "figure_study_architecture.png" + manuscript_pdf = manuscript / "figure_study_architecture.pdf" + canvas.save(manuscript_png, optimize=True) + canvas.save(manuscript_pdf, "PDF", resolution=300.0) + print(f"Wrote {out_png}") + print(f"Wrote {out_pdf}") + print(f"Wrote {manuscript_png}") + print(f"Wrote {manuscript_pdf}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/paper_runs/paper_04/scripts/generate_comparison_requests.py b/examples/paper_runs/paper_04/scripts/generate_comparison_requests.py index 30b19eb..2fa95cb 100644 --- a/examples/paper_runs/paper_04/scripts/generate_comparison_requests.py +++ b/examples/paper_runs/paper_04/scripts/generate_comparison_requests.py @@ -359,6 +359,48 @@ def _apply_measurement_noise(bits: list[int], p_meas: float, rng: random.Random) bits[i] ^= 1 +def _logical_x_indices(geom: SurfaceGeometry) -> list[int]: + y_mid = geom.distance // 2 + return [y_mid * (geom.distance - 1) + x for x in range(geom.distance - 1)] + + +def _error_indices(bits: list[int]) -> list[int]: + return [idx for idx, bit in enumerate(bits) if bit & 1] + + +def _parity_on(indices: list[int], bits: list[int]) -> int: + parity = 0 + for idx in indices: + if 0 <= idx < len(bits): + parity ^= bits[idx] & 1 + return parity + + +def _truth_record( + *, + code_id: str, + shot_index: int, + dataset_name: str, + geom: SurfaceGeometry, + x_bits: list[int], + z_bits: list[int], + truth_model: str, +) -> dict[str, Any]: + logical_indices = _logical_x_indices(geom) + return { + "code_id": code_id, + "round_index": shot_index, + "dataset": dataset_name, + "n_qubits": geom.n_data, + "logical_observable": "x_error_midline_parity", + "logical_indices": logical_indices, + "logical_truth": _parity_on(logical_indices, x_bits), + "x_error_indices": _error_indices(x_bits), + "z_error_indices": _error_indices(z_bits), + "truth_model": truth_model, + } + + def _simulate_repeated_round_events( *, geom: SurfaceGeometry, @@ -370,7 +412,7 @@ def _simulate_repeated_round_events( rng: random.Random, emit_x_events: bool, emit_z_events: bool, -) -> list[dict[str, Any]]: +) -> tuple[list[dict[str, Any]], list[int], list[int]]: x_bits = [0] * geom.n_data z_bits = [0] * geom.n_data prev_sx = [0] * geom.n_x @@ -395,7 +437,7 @@ def _simulate_repeated_round_events( prev_sx = sx prev_sz = sz - return events + return events, x_bits, z_bits def _digitize_periodic(value: float, period: float, width: float, bias: float = 0.0) -> int: @@ -493,7 +535,7 @@ def _simulate_gkp_repeated_round_events( rng: random.Random, emit_x_events: bool, emit_z_events: bool, -) -> list[dict[str, Any]]: +) -> tuple[list[dict[str, Any]], list[int], list[int]]: q_shift = [0.0] * geom.n_data p_shift = [0.0] * geom.n_data prev_sx = [0] * geom.n_x @@ -528,12 +570,16 @@ def _simulate_gkp_repeated_round_events( prev_sx = sx prev_sz = sz - return events + sqrt_pi = math.sqrt(math.pi) + x_bits = [_digitize_periodic(v, period=sqrt_pi, width=0.25 * sqrt_pi, bias=0.0) for v in q_shift] + z_bits = [_digitize_periodic(v, period=sqrt_pi, width=0.25 * sqrt_pi, bias=0.0) for v in p_shift] + return events, x_bits, z_bits def _write_dataset( path: Path, *, + truth_path: Path, dataset_name: str, code_id: str, shots: int, @@ -553,9 +599,9 @@ def _write_dataset( total_events = 0 rng = random.Random(seed) - with path.open("w", encoding="utf-8") as f: + with path.open("w", encoding="utf-8") as f, truth_path.open("w", encoding="utf-8") as tf: for shot_index in range(shots): - events = _simulate_repeated_round_events( + events, x_bits, z_bits = _simulate_repeated_round_events( geom=geom, rounds=rounds, p_gate=p_gate, @@ -596,10 +642,26 @@ def _write_dataset( }, } f.write(json.dumps(rec, separators=(",", ":")) + "\n") + tf.write( + json.dumps( + _truth_record( + code_id=code_id, + shot_index=shot_index, + dataset_name=dataset_name, + geom=geom, + x_bits=x_bits, + z_bits=z_bits, + truth_model="surface_final_pauli_state", + ), + separators=(",", ":"), + ) + + "\n" + ) return { "dataset": dataset_name, "request_file": path.name, + "truth_file": truth_path.name, "request_lines": shots, "avg_request_events": float(total_events) / float(max(shots, 1)), "nonempty_request_event_rate": float(nonempty) / float(max(shots, 1)), @@ -610,6 +672,7 @@ def _write_dataset( def _write_dataset_gkp( path: Path, *, + truth_path: Path, dataset_name: str, code_id: str, shots: int, @@ -630,9 +693,9 @@ def _write_dataset_gkp( total_events = 0 rng = random.Random(seed) - with path.open("w", encoding="utf-8") as f: + with path.open("w", encoding="utf-8") as f, truth_path.open("w", encoding="utf-8") as tf: for shot_index in range(shots): - events = _simulate_gkp_repeated_round_events( + events, x_bits, z_bits = _simulate_gkp_repeated_round_events( geom=geom, rounds=rounds, sigma_shift=sigma_shift, @@ -674,10 +737,26 @@ def _write_dataset_gkp( }, } f.write(json.dumps(rec, separators=(",", ":")) + "\n") + tf.write( + json.dumps( + _truth_record( + code_id=code_id, + shot_index=shot_index, + dataset_name=dataset_name, + geom=geom, + x_bits=x_bits, + z_bits=z_bits, + truth_model="gkp_reference_digitized_shift", + ), + separators=(",", ":"), + ) + + "\n" + ) return { "dataset": dataset_name, "request_file": path.name, + "truth_file": truth_path.name, "request_lines": shots, "avg_request_events": float(total_events) / float(max(shots, 1)), "nonempty_request_event_rate": float(nonempty) / float(max(shots, 1)), @@ -746,6 +825,7 @@ def main() -> int: rows.append( _write_dataset( out_dir / "decoder_requests_pennylane.ndjson", + truth_path=out_dir / "truth_pennylane.ndjson", dataset_name="pennylane", code_id=code_id, shots=args.shots, @@ -765,6 +845,7 @@ def main() -> int: rows.append( _write_dataset( out_dir / "decoder_requests_qiskit.ndjson", + truth_path=out_dir / "truth_qiskit.ndjson", dataset_name="qiskit", code_id=code_id, shots=args.shots, @@ -784,6 +865,7 @@ def main() -> int: rows.append( _write_dataset( out_dir / "decoder_requests_cirq.ndjson", + truth_path=out_dir / "truth_cirq.ndjson", dataset_name="cirq", code_id=code_id, shots=args.shots, @@ -803,6 +885,7 @@ def main() -> int: rows.append( _write_dataset( out_dir / "decoder_requests_lidmas_reference.ndjson", + truth_path=out_dir / "truth_lidmas_reference.ndjson", dataset_name="lidmas_reference", code_id=code_id, shots=args.shots, @@ -832,6 +915,7 @@ def main() -> int: rows.append( _write_dataset_gkp( out_dir / "decoder_requests_pennylane.ndjson", + truth_path=out_dir / "truth_pennylane.ndjson", dataset_name="pennylane", code_id=code_id, shots=args.shots, @@ -852,6 +936,7 @@ def main() -> int: rows.append( _write_dataset_gkp( out_dir / "decoder_requests_qiskit.ndjson", + truth_path=out_dir / "truth_qiskit.ndjson", dataset_name="qiskit", code_id=code_id, shots=args.shots, @@ -872,6 +957,7 @@ def main() -> int: rows.append( _write_dataset_gkp( out_dir / "decoder_requests_cirq.ndjson", + truth_path=out_dir / "truth_cirq.ndjson", dataset_name="cirq", code_id=code_id, shots=args.shots, @@ -892,6 +978,7 @@ def main() -> int: rows.append( _write_dataset_gkp( out_dir / "decoder_requests_lidmas_reference.ndjson", + truth_path=out_dir / "truth_lidmas_reference.ndjson", dataset_name="lidmas_reference", code_id=code_id, shots=args.shots, @@ -917,6 +1004,7 @@ def main() -> int: fieldnames=[ "dataset", "request_file", + "truth_file", "request_lines", "avg_request_events", "nonempty_request_event_rate", diff --git a/examples/paper_runs/paper_04/scripts/print_run_circuits.py b/examples/paper_runs/paper_04/scripts/print_run_circuits.py new file mode 100644 index 0000000..242918e --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/print_run_circuits.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 +"""Print circuit/logic snapshots used by paper_04 runs.""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path +from typing import Any + +from generate_comparison_requests import SurfaceGeometry, build_surface_geometry + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--run-root", + default="examples/paper_runs/paper_04/results/03_analysis/runs", + help="Root directory containing per-family runs.", + ) + parser.add_argument( + "--out-dir", + default="examples/paper_runs/paper_04/results/03_analysis/circuit_prints", + help="Output directory for printed circuit/logic text files.", + ) + return parser.parse_args() + + +def _read_json(path: Path) -> dict[str, Any]: + with path.open("r", encoding="utf-8") as f: + return json.load(f) + + +def _fmt_supports(geom: SurfaceGeometry) -> str: + lines: list[str] = [] + lines.append(f"distance={geom.distance}") + lines.append(f"n_data={geom.n_data}, n_x_checks={geom.n_x}, n_z_checks={geom.n_z}") + lines.append("") + lines.append("X-check supports (index: data-qubit list):") + for idx, support in enumerate(geom.x_supports): + lines.append(f" X{idx:02d}: {support}") + lines.append("") + lines.append("Z-check supports (index: data-qubit list):") + for idx, support in enumerate(geom.z_supports): + lines.append(f" Z{idx:02d}: {support}") + lines.append("") + return "\n".join(lines) + + +def _surface_qiskit_circuit_text(geom: SurfaceGeometry) -> str: + try: + from qiskit import QuantumCircuit # type: ignore + except Exception as exc: + return f"Qiskit unavailable in this environment: {exc}\n" + + n_total = geom.n_data + geom.n_x + geom.n_z + x_offset = geom.n_data + z_offset = geom.n_data + geom.n_x + qc = QuantumCircuit(n_total, name="surface_round") + + for c_idx, support in enumerate(geom.x_supports): + anc = x_offset + c_idx + qc.h(anc) + for dq in support: + qc.cx(anc, dq) + qc.h(anc) + + for c_idx, support in enumerate(geom.z_supports): + anc = z_offset + c_idx + for dq in support: + qc.cx(dq, anc) + + return str(qc.draw(output="text", fold=-1)) + + +def _surface_cirq_circuit_text(geom: SurfaceGeometry) -> str: + try: + import cirq # type: ignore + except Exception as exc: + return f"Cirq unavailable in this environment: {exc}\n" + + n_total = geom.n_data + geom.n_x + geom.n_z + x_offset = geom.n_data + z_offset = geom.n_data + geom.n_x + qubits = cirq.LineQubit.range(n_total) + x_anc = [qubits[x_offset + i] for i in range(geom.n_x)] + z_anc = [qubits[z_offset + i] for i in range(geom.n_z)] + + ops = [] + for c_idx, support in enumerate(geom.x_supports): + anc = x_anc[c_idx] + ops.append(cirq.H(anc)) + for dq in support: + ops.append(cirq.CNOT(anc, qubits[dq])) + ops.append(cirq.H(anc)) + + for c_idx, support in enumerate(geom.z_supports): + anc = z_anc[c_idx] + for dq in support: + ops.append(cirq.CNOT(qubits[dq], anc)) + + circuit = cirq.Circuit( + ops, + cirq.measure(*x_anc, key="mx"), + cirq.measure(*z_anc, key="mz"), + ) + return str(circuit) + + +def _surface_pennylane_circuit_text(geom: SurfaceGeometry) -> str: + try: + import pennylane as qml # type: ignore + except Exception as exc: + return f"PennyLane unavailable in this environment: {exc}\n" + + n_total = geom.n_data + geom.n_x + geom.n_z + x_offset = geom.n_data + z_offset = geom.n_data + geom.n_x + dev = qml.device("default.clifford", wires=n_total) + + @qml.set_shots(shots=1) + @qml.qnode(dev) + def circuit(): + for c_idx, support in enumerate(geom.x_supports): + anc = x_offset + c_idx + qml.Hadamard(anc) + for dq in support: + qml.CNOT(wires=[anc, dq]) + qml.Hadamard(anc) + + for c_idx, support in enumerate(geom.z_supports): + anc = z_offset + c_idx + for dq in support: + qml.CNOT(wires=[dq, anc]) + + measures = [qml.sample(qml.PauliZ(x_offset + i)) for i in range(geom.n_x)] + measures.extend(qml.sample(qml.PauliZ(z_offset + i)) for i in range(geom.n_z)) + return measures + + try: + drawer = qml.draw(circuit, expansion_strategy="device") + except TypeError: + drawer = qml.draw(circuit) + return drawer() + + +def _gkp_logic_text(variant: str, geom: SurfaceGeometry) -> str: + lines: list[str] = [] + lines.append(f"GKP digitized logic variant: {variant}") + lines.append("") + lines.append("No framework-native gate circuit is built for GKP in this run.") + lines.append("Instead, repeated-round q/p shift states are projected through check supports") + lines.append("and digitized into X/Z syndrome bits with variant-specific rules.") + lines.append("") + lines.append("Digitization rules:") + lines.append(" - pennylane: periodic threshold with small bias") + lines.append(" - qiskit: rounded scaled value with Gaussian perturbation") + lines.append(" - cirq: sinusoidal phase-sign rule") + lines.append(" - lidmas_reference: periodic threshold without framework bias") + lines.append("") + lines.append(_fmt_supports(geom)) + return "\n".join(lines) + + +def _run_meta_text(summary: dict[str, Any]) -> str: + keys = [ + "shots", + "rounds", + "distance", + "code_family", + "n_qubits", + "n_x_checks", + "n_z_checks", + "error_rate", + "sigma", + "seed", + "emit_x_events", + "emit_z_events", + "pennylane_enabled", + "qiskit_enabled", + "cirq_enabled", + ] + lines = ["Run metadata:"] + for key in keys: + lines.append(f" {key}: {summary.get(key)}") + return "\n".join(lines) + "\n" + + +def main() -> int: + args = parse_args() + run_root = Path(args.run_root) + out_dir = Path(args.out_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + # Use current manuscript run geometry (distance from summary). + surf_summary = _read_json(run_root / "surface/01_generate_comparison_requests/summary_generation.json") + gkp_summary = _read_json(run_root / "gkp/01_generate_comparison_requests/summary_generation.json") + + surf_geom = build_surface_geometry(int(surf_summary["distance"])) + gkp_geom = build_surface_geometry(int(gkp_summary["distance"])) + + # Shared metadata files. + (out_dir / "surface_run_metadata.txt").write_text(_run_meta_text(surf_summary), encoding="utf-8") + (out_dir / "gkp_run_metadata.txt").write_text(_run_meta_text(gkp_summary), encoding="utf-8") + (out_dir / "surface_check_supports.txt").write_text(_fmt_supports(surf_geom), encoding="utf-8") + (out_dir / "gkp_check_supports.txt").write_text(_fmt_supports(gkp_geom), encoding="utf-8") + + # Surface family circuits. + (out_dir / "surface_pennylane_circuit.txt").write_text( + _surface_pennylane_circuit_text(surf_geom), encoding="utf-8" + ) + (out_dir / "surface_qiskit_circuit.txt").write_text( + _surface_qiskit_circuit_text(surf_geom), encoding="utf-8" + ) + (out_dir / "surface_cirq_circuit.txt").write_text( + _surface_cirq_circuit_text(surf_geom), encoding="utf-8" + ) + (out_dir / "surface_lidmas_reference_logic.txt").write_text( + "LiDMaS+ reference path uses the classical parity-check sampler.\n\n" + _fmt_supports(surf_geom), + encoding="utf-8", + ) + + # GKP family logic (digitized, not gate-circuit objects). + for variant in ("pennylane", "qiskit", "cirq", "lidmas_reference"): + (out_dir / f"gkp_{variant}_digitized_logic.txt").write_text( + _gkp_logic_text(variant, gkp_geom), encoding="utf-8" + ) + + print(f"Wrote circuit/logic printouts to {out_dir}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/paper_runs/paper_04/scripts/render_2d_syndrome_layouts.py b/examples/paper_runs/paper_04/scripts/render_2d_syndrome_layouts.py new file mode 100644 index 0000000..9eb7764 --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/render_2d_syndrome_layouts.py @@ -0,0 +1,501 @@ +#!/usr/bin/env python3 +"""Render 2D syndrome-layout figures for paper_04 runs.""" + +from __future__ import annotations + +import argparse +import json +import math +from pathlib import Path +from typing import Any + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +from matplotlib.patches import Circle, Patch, Polygon + +from generate_comparison_requests import SurfaceGeometry, build_surface_geometry + +STYLE = { + "bg": "#FFFFFF", + "support": "#111827", + "data": "#DC2626", + "data_text": "#7F1D1D", + "x": "#2B7BBB", + "x_text": "#1D4ED8", + "x_stab": "#AFC7F2", + "z": "#4DA64D", + "z_text": "#166534", + "z_stab": "#B9DDB4", +} + +LABEL_BBOX = {"facecolor": "white", "edgecolor": "none", "alpha": 0.74, "pad": 0.04} + +plt.rcParams.update( + { + "font.family": "DejaVu Sans", + "axes.facecolor": STYLE["bg"], + "figure.facecolor": STYLE["bg"], + "savefig.facecolor": STYLE["bg"], + } +) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--run-root", + default="examples/paper_runs/paper_04/results/03_analysis/runs", + help="Root directory containing per-family run outputs.", + ) + parser.add_argument( + "--out-dir", + default="examples/paper_runs/paper_04/results/03_analysis/syndrome_layout_figures", + help="Output directory for 2D syndrome layout figures.", + ) + return parser.parse_args() + + +def _load_summary(path: Path) -> dict[str, Any]: + with path.open("r", encoding="utf-8") as f: + return json.load(f) + + +def _layout_coords(geom: SurfaceGeometry) -> tuple[dict[int, tuple[float, float]], dict[int, tuple[float, float]], dict[int, tuple[float, float]]]: + d = geom.distance + data: dict[int, tuple[float, float]] = {} + x_checks: dict[int, tuple[float, float]] = {} + z_checks: dict[int, tuple[float, float]] = {} + + idx = 0 + # Horizontal data qubits h(x, y): between X(x, y) and X(x+1, y) + for y in range(d): + for x in range(d - 1): + data[idx] = (x + 0.5, float(y)) + idx += 1 + # Vertical data qubits v(x, y): between X(x, y) and X(x, y+1) + for y in range(d - 1): + for x in range(d): + data[idx] = (float(x), y + 0.5) + idx += 1 + + for y in range(d): + for x in range(d): + x_idx = y * d + x + x_checks[x_idx] = (float(x), float(y)) + + z_idx = 0 + for y in range(d - 1): + for x in range(d - 1): + z_checks[z_idx] = (x + 0.5, y + 0.5) + z_idx += 1 + + return data, x_checks, z_checks + + +def _ordered_support_points(points: list[tuple[float, float]]) -> list[tuple[float, float]]: + if len(points) <= 2: + return points + cx = sum(x for x, _ in points) / len(points) + cy = sum(y for _, y in points) / len(points) + return sorted(points, key=lambda p: math.atan2(p[1] - cy, p[0] - cx)) + + +def _draw_stabilizer_patch(ax: Any, points: list[tuple[float, float]], color: str, zorder: int) -> None: + ordered = _ordered_support_points(points) + if len(ordered) >= 3: + ax.add_patch( + Polygon( + ordered, + closed=True, + facecolor=color, + edgecolor="none", + alpha=0.55, + zorder=zorder, + ) + ) + elif len(ordered) == 2: + (x0, y0), (x1, y1) = ordered + ax.plot([x0, x1], [y0, y1], color=color, linewidth=12.0, alpha=0.40, solid_capstyle="round", zorder=zorder) + + +def _draw_panel( + ax: Any, + geom: SurfaceGeometry, + *, + show_x: bool, + show_z: bool, + indexed: bool, + title: str, +) -> None: + data_coords, x_coords, z_coords = _layout_coords(geom) + d = geom.distance + + def _is_int(v: float) -> bool: + return abs(v - round(v)) < 1e-9 + + # Draw translucent stabilizer regions first, then overlay dashed supports. + if show_x: + for support in geom.x_supports: + _draw_stabilizer_patch(ax, [data_coords[dq] for dq in support], STYLE["x_stab"], zorder=0) + + if show_z: + for support in geom.z_supports: + _draw_stabilizer_patch(ax, [data_coords[dq] for dq in support], STYLE["z_stab"], zorder=0) + + if show_x: + for x_idx, support in enumerate(geom.x_supports): + cx, cy = x_coords[x_idx] + for dq in support: + qx, qy = data_coords[dq] + ax.plot( + [cx, qx], + [cy, qy], + color=STYLE["support"], + linewidth=0.85, + alpha=0.72, + linestyle=(0, (2.5, 2.5)), + zorder=1, + ) + + if show_z: + for z_idx, support in enumerate(geom.z_supports): + cx, cy = z_coords[z_idx] + for dq in support: + qx, qy = data_coords[dq] + ax.plot( + [cx, qx], + [cy, qy], + color=STYLE["support"], + linewidth=0.85, + alpha=0.72, + linestyle=(0, (2.5, 2.5)), + zorder=1, + ) + + # Data qubits. + for dq, (qx, qy) in data_coords.items(): + ax.scatter(qx, qy, s=54, color=STYLE["data"], edgecolors="white", linewidths=0.8, zorder=5) + if indexed: + # Route data labels by edge orientation to reduce overlap with check labels. + if _is_int(qy): # horizontal data edge + ddx, ddy = (0.06, -0.08) + else: # vertical data edge + ddx, ddy = (0.07, 0.08) + ax.text( + qx + ddx, + qy + ddy, + f"D{dq}", + fontsize=6.0, + color=STYLE["data_text"], + zorder=6, + bbox=LABEL_BBOX, + ) + + # X checks. + if show_x: + xs = [x for x, _ in x_coords.values()] + ys = [y for _, y in x_coords.values()] + ax.scatter( + xs, + ys, + s=72, + marker="o", + color=STYLE["x"], + edgecolors="white", + linewidths=0.8, + zorder=5, + label="X checks", + ) + if indexed: + for x_idx, (cx, cy) in x_coords.items(): + # Keep X labels above checks and nudge outer columns inward. + xdx = -0.13 + if cx <= 0.0: + xdx = -0.10 + elif cx >= d - 1: + xdx = -0.16 + ax.text( + cx + xdx, + cy - 0.19, + f"X{x_idx:02d}", + fontsize=6.0, + color=STYLE["x_text"], + zorder=5, + bbox=LABEL_BBOX, + ) + + # Z checks. + if show_z: + xs = [x for x, _ in z_coords.values()] + ys = [y for _, y in z_coords.values()] + ax.scatter( + xs, + ys, + s=72, + marker="o", + color=STYLE["z"], + edgecolors="white", + linewidths=0.8, + zorder=5, + label="Z checks", + ) + if indexed: + for z_idx, (cx, cy) in z_coords.items(): + # Place Z labels below checks, alternating side to avoid dense runs. + zdx = -0.15 if (z_idx % 2 == 0) else 0.05 + if cx <= 0.7: + zdx = 0.04 + elif cx >= d - 1.3: + zdx = -0.18 + ax.text( + cx + zdx, + cy + 0.13, + f"Z{z_idx:02d}", + fontsize=6.0, + color=STYLE["z_text"], + zorder=5, + bbox=LABEL_BBOX, + ) + ax.set_xlim(-0.60, d - 1 + 0.60) + ax.set_ylim(d - 1 + 0.60, -0.60) + ax.set_aspect("equal") + ax.set_xticks([]) + ax.set_yticks([]) + ax.grid(False) + for spine in ax.spines.values(): + spine.set_visible(False) + ax.set_title(title, fontsize=12, pad=10) + ax.set_xlabel("") + ax.set_ylabel("") + + +def _save_all(fig: Any, out_base: Path) -> None: + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_base.with_suffix(ext), bbox_inches="tight", dpi=360) + + +def render_surface_figures(geom: SurfaceGeometry, out_dir: Path) -> None: + # Compact manuscript-facing view. + fig, axes = plt.subplots(1, 3, figsize=(15.8, 5.3), constrained_layout=True) + _draw_panel(axes[0], geom, show_x=True, show_z=False, indexed=False, title="Surface: X-check supports") + _draw_panel(axes[1], geom, show_x=False, show_z=True, indexed=False, title="Surface: Z-check supports") + _draw_panel(axes[2], geom, show_x=True, show_z=True, indexed=False, title="Surface: Combined support graph") + + handles = [ + Line2D([0], [0], marker="o", color="none", markerfacecolor=STYLE["data"], markeredgecolor="white", label="Data qubit"), + Line2D([0], [0], marker="o", color="none", markerfacecolor=STYLE["x"], markeredgecolor="white", label="X ancilla qubit"), + Line2D([0], [0], marker="o", color="none", markerfacecolor=STYLE["z"], markeredgecolor="white", label="Z ancilla qubit"), + Patch(facecolor=STYLE["x_stab"], edgecolor="none", alpha=0.55, label="X stabilizer"), + Patch(facecolor=STYLE["z_stab"], edgecolor="none", alpha=0.55, label="Z stabilizer"), + ] + fig.legend(handles=handles, loc="upper center", ncol=5, frameon=False, bbox_to_anchor=(0.5, 1.13), fontsize=10) + _save_all(fig, out_dir / "figure_surface_2d_syndrome_layout") + plt.close(fig) + + # Indexed engineering/debug view. + fig_idx, ax_idx = plt.subplots(figsize=(7.8, 6.8), constrained_layout=True) + _draw_panel( + ax_idx, + geom, + show_x=True, + show_z=True, + indexed=True, + title="Surface: Indexed support graph (data/check ids)", + ) + _save_all(fig_idx, out_dir / "figure_surface_2d_syndrome_layout_indexed") + plt.close(fig_idx) + + # Indexed manuscript-facing triptych view (X-only, Z-only, combined). + fig_all_idx, axes_all_idx = plt.subplots(1, 3, figsize=(17.2, 5.7), constrained_layout=True) + _draw_panel( + axes_all_idx[0], + geom, + show_x=True, + show_z=False, + indexed=True, + title="Surface: X-check supports (indexed)", + ) + _draw_panel( + axes_all_idx[1], + geom, + show_x=False, + show_z=True, + indexed=True, + title="Surface: Z-check supports (indexed)", + ) + _draw_panel( + axes_all_idx[2], + geom, + show_x=True, + show_z=True, + indexed=True, + title="Surface: Combined support graph (indexed)", + ) + fig_all_idx.legend(handles=handles, loc="upper center", ncol=5, frameon=False, bbox_to_anchor=(0.5, 1.13), fontsize=10) + _save_all(fig_all_idx, out_dir / "figure_surface_2d_syndrome_layout_all_indexed") + plt.close(fig_all_idx) + + +def render_gkp_figures(geom: SurfaceGeometry, out_dir: Path) -> None: + # Same topology as surface, with GKP annotation. + fig, ax = plt.subplots(figsize=(7.6, 6.4), constrained_layout=True) + _draw_panel( + ax, + geom, + show_x=True, + show_z=True, + indexed=False, + title=( + "Digitized-GKP outer support topology\n" + "(shared by PennyLane / Qiskit / Cirq / LiDMaS+ variants)" + ), + ) + _save_all(fig, out_dir / "figure_gkp_2d_syndrome_layout") + plt.close(fig) + + # Indexed engineering/debug view. + fig_idx, ax_idx = plt.subplots(figsize=(8.4, 7.2), constrained_layout=True) + _draw_panel( + ax_idx, + geom, + show_x=True, + show_z=True, + indexed=True, + title=( + "Digitized-GKP outer support topology (indexed)\n" + "(shared by PennyLane / Qiskit / Cirq / LiDMaS+ variants)" + ), + ) + _save_all(fig_idx, out_dir / "figure_gkp_2d_syndrome_layout_indexed") + plt.close(fig_idx) + + +def _draw_gkp_phase_space(ax: Any, *, indexed: bool = False) -> None: + sqrt_pi = math.sqrt(math.pi) + lim = 2.5 * sqrt_pi + + # Stabilizer decision boundaries. + for k in range(-2, 3): + x = (k + 0.5) * sqrt_pi + y = (k + 0.5) * sqrt_pi + ax.axvline(x, color="#94A3B8", linewidth=1.0, linestyle="--", alpha=0.8, zorder=1) + ax.axhline(y, color="#94A3B8", linewidth=1.0, linestyle="--", alpha=0.8, zorder=1) + + # GKP peak lattice (conceptual envelope). + for i in range(-2, 3): + for j in range(-2, 3): + q = i * sqrt_pi + p = j * sqrt_pi + w = 1.0 if (i == 0 and j == 0) else 0.6 + ax.scatter(q, p, s=36 * w, color="#2563EB", alpha=0.78, zorder=3) + ax.add_patch(Circle((q, p), radius=0.10 * sqrt_pi, edgecolor="#2563EB", facecolor="none", alpha=0.35, zorder=2)) + if indexed: + ax.text( + q + 0.08 * sqrt_pi, + p + 0.08 * sqrt_pi, + f"({i},{j})", + fontsize=6.3, + color="#1D4ED8", + zorder=4, + ) + + ax.set_xlim(-lim, lim) + ax.set_ylim(-lim, lim) + ax.set_aspect("equal") + ax.grid(color="#CBD5E1", alpha=0.45, linewidth=0.6) + ax.set_xlabel("q quadrature") + ax.set_ylabel("p quadrature") + title = "Inner GKP code (indexed phase-space cell)" if indexed else "Inner GKP code (single-mode phase-space cell)" + ax.set_title( + title + "\n" + r"Peak spacing $\sqrt{\pi}$; dashed lines indicate digitization boundaries.", + fontsize=10, + pad=10, + ) + + +def _draw_outer_minimap(ax: Any, geom: SurfaceGeometry, *, indexed: bool = False) -> None: + _draw_panel( + ax, + geom, + show_x=True, + show_z=True, + indexed=indexed, + title="Outer code support graph (indexed)" if indexed else "Outer code support graph", + ) + ax.set_xlabel("") + ax.set_ylabel("") + + +def render_gkp_inner_outer_figures(geom: SurfaceGeometry, out_dir: Path) -> None: + # Standalone inner-code conceptual figure. + fig_inner, ax_inner = plt.subplots(figsize=(6.2, 5.5), constrained_layout=True) + _draw_gkp_phase_space(ax_inner) + _save_all(fig_inner, out_dir / "figure_gkp_inner_code_phase_space") + plt.close(fig_inner) + + fig_inner_idx, ax_inner_idx = plt.subplots(figsize=(6.4, 5.7), constrained_layout=True) + _draw_gkp_phase_space(ax_inner_idx, indexed=True) + _save_all(fig_inner_idx, out_dir / "figure_gkp_inner_code_phase_space_indexed") + plt.close(fig_inner_idx) + + # Concatenated schematic: inner GKP -> outer support graph. + fig_cat, axes = plt.subplots(1, 2, figsize=(12.4, 5.4), constrained_layout=True) + _draw_gkp_phase_space(axes[0]) + _draw_outer_minimap(axes[1], geom) + axes[0].set_title("Inner code: GKP digitization space", fontsize=10, pad=10) + axes[1].set_title("Outer code: support topology", fontsize=10, pad=10) + + # Cross-panel mapping annotation. + axes[0].annotate( + "", + xy=(1.03, 0.5), + xycoords="axes fraction", + xytext=(-0.03, 0.5), + textcoords=axes[1].transAxes, + arrowprops={"arrowstyle": "->", "color": "#0F172A", "lw": 1.6}, + ) + _save_all(fig_cat, out_dir / "figure_gkp_concatenated_inner_outer") + plt.close(fig_cat) + + # Indexed concatenated schematic (right panel indexed). + fig_cat_idx, axes_idx = plt.subplots(1, 2, figsize=(13.2, 6.2), constrained_layout=True) + _draw_gkp_phase_space(axes_idx[0], indexed=True) + _draw_outer_minimap(axes_idx[1], geom, indexed=True) + axes_idx[0].set_title("Inner code: GKP digitization space", fontsize=10, pad=10) + axes_idx[1].set_title("Outer code: support topology (indexed)", fontsize=10, pad=10) + + axes_idx[0].annotate( + "", + xy=(1.03, 0.5), + xycoords="axes fraction", + xytext=(-0.03, 0.5), + textcoords=axes_idx[1].transAxes, + arrowprops={"arrowstyle": "->", "color": "#0F172A", "lw": 1.6}, + ) + _save_all(fig_cat_idx, out_dir / "figure_gkp_concatenated_inner_outer_indexed") + plt.close(fig_cat_idx) + + +def main() -> int: + args = parse_args() + run_root = Path(args.run_root) + out_dir = Path(args.out_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + surface_summary = _load_summary(run_root / "surface/01_generate_comparison_requests/summary_generation.json") + gkp_summary = _load_summary(run_root / "gkp/01_generate_comparison_requests/summary_generation.json") + + surface_geom = build_surface_geometry(int(surface_summary["distance"])) + gkp_geom = build_surface_geometry(int(gkp_summary["distance"])) + + render_surface_figures(surface_geom, out_dir) + render_gkp_figures(gkp_geom, out_dir) + render_gkp_inner_outer_figures(gkp_geom, out_dir) + print(f"Wrote 2D syndrome layout figures to {out_dir}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/paper_runs/paper_04/scripts/render_3d_lattice_preview.py b/examples/paper_runs/paper_04/scripts/render_3d_lattice_preview.py new file mode 100644 index 0000000..bb354c1 --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/render_3d_lattice_preview.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python3 +"""Render a 3D preview of the indexed outer-code lattice for paper_04.""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path +from typing import Any + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +from matplotlib.patches import Patch +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + +from generate_comparison_requests import SurfaceGeometry, build_surface_geometry +from render_2d_syndrome_layouts import STYLE, _layout_coords, _ordered_support_points + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--summary", + default="examples/paper_runs/paper_04/results/03_analysis/runs/surface/01_generate_comparison_requests/summary_generation.json", + help="Generation summary JSON used to recover the run distance.", + ) + parser.add_argument( + "--out-dir", + default="examples/paper_runs/paper_04/results/03_analysis/syndrome_layout_figures", + help="Output directory for the 3D preview figure.", + ) + return parser.parse_args() + + +def _load_summary(path: Path) -> dict[str, Any]: + with path.open("r", encoding="utf-8") as f: + return json.load(f) + + +def _patch_3d(points: list[tuple[float, float]], z: float) -> list[tuple[float, float, float]]: + return [(x, y, z) for x, y in _ordered_support_points(points)] + + +def render_3d_lattice(geom: SurfaceGeometry, out_dir: Path) -> None: + data_coords, x_coords, z_coords = _layout_coords(geom) + fig = plt.figure(figsize=(9.2, 7.2), dpi=320) + ax = fig.add_subplot(111, projection="3d") + + z_data = 0.0 + z_x = 0.42 + z_z = -0.42 + z_x_patch = 0.18 + z_z_patch = -0.18 + + # Stabilizer sheets. + x_polys = [_patch_3d([data_coords[dq] for dq in support], z_x_patch) for support in geom.x_supports] + z_polys = [_patch_3d([data_coords[dq] for dq in support], z_z_patch) for support in geom.z_supports] + ax.add_collection3d( + Poly3DCollection(x_polys, facecolors=STYLE["x_stab"], edgecolors="none", alpha=0.42, zorder=0) + ) + ax.add_collection3d( + Poly3DCollection(z_polys, facecolors=STYLE["z_stab"], edgecolors="none", alpha=0.46, zorder=0) + ) + + # Dashed support links from ancillas to data qubits. + for x_idx, support in enumerate(geom.x_supports): + cx, cy = x_coords[x_idx] + for dq in support: + qx, qy = data_coords[dq] + ax.plot( + [cx, qx], + [cy, qy], + [z_x, z_data], + color=STYLE["support"], + linewidth=0.8, + alpha=0.65, + linestyle=(0, (2.5, 2.5)), + ) + + for z_idx, support in enumerate(geom.z_supports): + cx, cy = z_coords[z_idx] + for dq in support: + qx, qy = data_coords[dq] + ax.plot( + [cx, qx], + [cy, qy], + [z_z, z_data], + color=STYLE["support"], + linewidth=0.8, + alpha=0.65, + linestyle=(0, (2.5, 2.5)), + ) + + # Nodes. + dx = [x for x, _ in data_coords.values()] + dy = [y for _, y in data_coords.values()] + ax.scatter(dx, dy, [z_data] * len(dx), s=42, color=STYLE["data"], edgecolors="white", linewidths=0.7, depthshade=False) + + xx = [x for x, _ in x_coords.values()] + xy = [y for _, y in x_coords.values()] + ax.scatter(xx, xy, [z_x] * len(xx), s=54, color=STYLE["x"], edgecolors="white", linewidths=0.7, depthshade=False) + + zx = [x for x, _ in z_coords.values()] + zy = [y for _, y in z_coords.values()] + ax.scatter(zx, zy, [z_z] * len(zx), s=54, color=STYLE["z"], edgecolors="white", linewidths=0.7, depthshade=False) + + # Sparse labels keep the preview readable in perspective. + for dq, (x, y) in data_coords.items(): + ax.text(x + 0.03, y + 0.03, z_data + 0.03, f"D{dq}", fontsize=5.0, color=STYLE["data_text"]) + for x_idx, (x, y) in x_coords.items(): + ax.text(x - 0.08, y - 0.06, z_x + 0.04, f"X{x_idx:02d}", fontsize=5.2, color=STYLE["x_text"]) + for z_idx, (x, y) in z_coords.items(): + ax.text(x + 0.03, y + 0.03, z_z - 0.03, f"Z{z_idx:02d}", fontsize=5.2, color=STYLE["z_text"]) + + d = geom.distance + ax.set_xlim(-0.55, d - 1 + 0.55) + ax.set_ylim(d - 1 + 0.55, -0.55) + ax.set_zlim(-0.72, 0.72) + ax.set_box_aspect((1.0, 1.0, 0.34)) + ax.view_init(elev=29, azim=-55) + ax.set_xlabel("lattice x", labelpad=8) + ax.set_ylabel("lattice y", labelpad=8) + ax.set_zlabel("syndrome layer", labelpad=8) + ax.set_zticks([z_z, z_data, z_x]) + ax.set_zticklabels(["Z", "data", "X"]) + ax.set_title("3D indexed outer-code lattice preview", pad=16) + ax.grid(alpha=0.16) + + legend_handles = [ + Line2D([0], [0], marker="o", color="none", markerfacecolor=STYLE["data"], markeredgecolor="white", label="Data qubit"), + Line2D([0], [0], marker="o", color="none", markerfacecolor=STYLE["x"], markeredgecolor="white", label="X ancilla"), + Line2D([0], [0], marker="o", color="none", markerfacecolor=STYLE["z"], markeredgecolor="white", label="Z ancilla"), + Patch(facecolor=STYLE["x_stab"], edgecolor="none", alpha=0.42, label="X stabilizer sheet"), + Patch(facecolor=STYLE["z_stab"], edgecolor="none", alpha=0.46, label="Z stabilizer sheet"), + ] + fig.legend(handles=legend_handles, loc="upper center", ncol=5, frameon=False, bbox_to_anchor=(0.5, 0.98), fontsize=9) + fig.subplots_adjust(top=0.88, left=0.02, right=0.98, bottom=0.02) + + out_dir.mkdir(parents=True, exist_ok=True) + out_base = out_dir / "figure_3d_lattice_preview" + for ext in (".png", ".pdf", ".svg"): + fig.savefig(out_base.with_suffix(ext), bbox_inches="tight", dpi=320) + plt.close(fig) + + +def main() -> int: + args = parse_args() + summary = _load_summary(Path(args.summary)) + geom = build_surface_geometry(int(summary["distance"])) + render_3d_lattice(geom, Path(args.out_dir)) + print(f"Wrote 3D lattice preview to {args.out_dir}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/paper_runs/paper_04/scripts/render_circuit_atlases.py b/examples/paper_runs/paper_04/scripts/render_circuit_atlases.py new file mode 100644 index 0000000..d6f380f --- /dev/null +++ b/examples/paper_runs/paper_04/scripts/render_circuit_atlases.py @@ -0,0 +1,199 @@ +#!/usr/bin/env python3 +"""Assemble circuit-only atlas figures for paper_04.""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +from PIL import Image, ImageDraw, ImageFont + + +PRINT_ROOT = Path("examples/paper_runs/paper_04/results/03_analysis/circuit_prints") +OUT_ROOT = Path("examples/paper_runs/paper_04/results/03_analysis/circuit_figures") + +TEXT_FONT = Path("./.venv/lib/python3.14/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSansMono.ttf") +BOLD_FONT = Path("./.venv/lib/python3.14/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSansMono-Bold.ttf") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--print-dir", + default=str(PRINT_ROOT), + help="Directory containing circuit text snapshots.", + ) + parser.add_argument( + "--out-dir", + default=str(OUT_ROOT), + help="Directory for grouped atlas outputs.", + ) + return parser.parse_args() + + +def _font(size: int, *, bold: bool = False) -> ImageFont.FreeTypeFont | ImageFont.ImageFont: + candidates = [ + BOLD_FONT if bold else TEXT_FONT, + Path("/System/Library/Fonts/Menlo.ttc"), + Path("/System/Library/Fonts/SFNSMono.ttf"), + Path("/System/Library/Fonts/Supplemental/Andale Mono.ttf"), + ] + for path in candidates: + if path.exists(): + try: + return ImageFont.truetype(str(path), size=size) + except OSError: + continue + return ImageFont.load_default() + + +def _read_blocks(path: Path) -> list[list[str]]: + text = path.read_text(encoding="utf-8").rstrip() + return [block.splitlines() for block in text.split("\n\n") if block.strip()] + + +def _text_size(lines: list[str], font: ImageFont.ImageFont, line_gap: int) -> tuple[int, int, int]: + sample = "M" + bbox = font.getbbox(sample) + line_h = bbox[3] - bbox[1] + line_gap + max_w = 1 + for line in lines: + line_bbox = font.getbbox(line or " ") + max_w = max(max_w, line_bbox[2] - line_bbox[0]) + return max_w, line_h * len(lines), line_h + + +def _render_text_panel( + lines: list[str], + *, + label: str, + font_size: int, + pad: int = 24, + line_gap: int = 3, +) -> Image.Image: + font = _font(font_size) + label_font = _font(max(13, font_size + 1), bold=True) + text_w, text_h, line_h = _text_size(lines, font, line_gap) + label_h = label_font.getbbox(label)[3] - label_font.getbbox(label)[1] + 16 + image = Image.new("RGB", (text_w + pad * 2, text_h + label_h + pad * 2), "#FFFFFF") + draw = ImageDraw.Draw(image) + draw.rounded_rectangle( + (0, 0, image.width - 1, image.height - 1), + radius=10, + fill="#FFFFFF", + outline="#CBD5E1", + width=2, + ) + draw.text((pad, pad - 2), label, fill="#111827", font=label_font) + y = pad + label_h + for line in lines: + draw.text((pad, y), line, fill="#111827", font=font) + y += line_h + return image + + +def _paste(canvas: Image.Image, panel: Image.Image, x: int, y: int) -> None: + canvas.paste(panel, (x, y)) + + +def _widen_to(panel: Image.Image, target_width: int) -> Image.Image: + if panel.width == target_width: + return panel + return panel.resize((target_width, panel.height), Image.Resampling.BICUBIC) + + +def _pennylane_row(print_dir: Path) -> tuple[str, list[Image.Image]]: + blocks = _read_blocks(print_dir / "surface_pennylane_circuit.txt") + panels = [ + _render_text_panel( + block, + label=f"PennyLane block {idx}/6, wires 0-80", + font_size=13, + pad=20, + line_gap=3, + ) + for idx, block in enumerate(blocks, start=1) + ] + return "PennyLane surface circuit unfolded left-to-right", panels + + +def _single_block_row(print_dir: Path, filename: str, title: str, font_size: int) -> tuple[str, list[Image.Image]]: + block = _read_blocks(print_dir / filename)[0] + panel = _render_text_panel(block, label=title, font_size=font_size, pad=22, line_gap=2) + return title, [panel] + + +def _compose_surface_atlas(print_dir: Path, out_dir: Path) -> None: + rows = [ + _pennylane_row(print_dir), + _single_block_row(print_dir, "surface_qiskit_circuit.txt", "Qiskit surface circuit", 8), + _single_block_row(print_dir, "surface_cirq_circuit.txt", "Cirq surface circuit", 10), + ] + + margin = 70 + gap = 30 + row_gap = 72 + title_h = 86 + note_h = 58 + title_font = _font(34, bold=True) + row_font = _font(22, bold=True) + note_font = _font(18) + + penny_title, penny_panels = rows[0] + penny_w = sum(panel.width for panel in penny_panels) + gap * (len(penny_panels) - 1) + rows = [rows[0]] + [(title, [_widen_to(panels[0], penny_w)]) for title, panels in rows[1:]] + max_w = penny_w + + width = margin * 2 + max_w + height = margin + title_h + for _, panels in rows: + height += row_font.getbbox("M")[3] - row_font.getbbox("M")[1] + 22 + height += max(panel.height for panel in panels) + height += row_gap + height += note_h + margin - row_gap + + canvas = Image.new("RGB", (width, height), "#FFFFFF") + draw = ImageDraw.Draw(canvas) + y = margin + draw.text((margin, y), "Circuit-only atlas for surface-family runs", fill="#111827", font=title_font) + y += title_h + + for row_title, panels in rows: + draw.text((margin, y), row_title, fill="#1F2937", font=row_font) + y += row_font.getbbox("M")[3] - row_font.getbbox("M")[1] + 22 + + if len(panels) == 1: + x = margin + _paste(canvas, panels[0], x, y) + y += panels[0].height + row_gap + continue + + x = margin + for panel in panels: + _paste(canvas, panel, x, y) + x += panel.width + gap + y += max(panel.height for panel in panels) + row_gap + + note = ( + "Only framework-rendered circuit drawings are shown here; support maps, metadata, " + "and GKP digitization tables are reported separately in the appendix tables." + ) + draw.text((margin, height - margin - note_h + 12), note, fill="#374151", font=note_font) + + out_dir.mkdir(parents=True, exist_ok=True) + png = out_dir / "figure_surface_circuit_atlas.png" + pdf = out_dir / "figure_surface_circuit_atlas.pdf" + canvas.save(png, optimize=True) + canvas.save(pdf, "PDF", resolution=300.0) + print(f"Wrote {png}") + print(f"Wrote {pdf}") + + +def main() -> int: + args = parse_args() + _compose_surface_atlas(Path(args.print_dir), Path(args.out_dir)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main())