From 3cbefd72a1c444ba9195cf7c22f8715c89132960 Mon Sep 17 00:00:00 2001 From: rosspeili Date: Tue, 9 Jun 2026 14:34:02 +0300 Subject: [PATCH 1/2] Add global rotations tutorial reproducing ideal Fig. 4a (issue #119). Introduces a docs demo with hypercube [[7,1,3]] flow verification, [[15,1,3]] stabilizer prep, and smoke tests for qualitative curve shapes. --- .gitignore | 3 + docs/demos/global_rotations.ipynb | 255 ++++++++++++++ docs/demos/utils/global_rotation.py | 420 ++++++++++++++++++++++++ mkdocs.yml | 1 + test/unit/demos/test_global_rotation.py | 45 +++ 5 files changed, 724 insertions(+) create mode 100644 docs/demos/global_rotations.ipynb create mode 100644 docs/demos/utils/global_rotation.py create mode 100644 test/unit/demos/test_global_rotation.py diff --git a/.gitignore b/.gitignore index 36f7ed60..98912e7c 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,9 @@ .vscode experiments/ +# Local reference material (papers, notes) — not part of the repo +*.pdf + # mkdocs site/ docs/reference/ diff --git a/docs/demos/global_rotations.ipynb b/docs/demos/global_rotations.ipynb new file mode 100644 index 00000000..dbee7a8b --- /dev/null +++ b/docs/demos/global_rotations.ipynb @@ -0,0 +1,255 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "71569fb7", + "metadata": {}, + "source": [ + "# Global rotations on colour codes\n", + "\n", + "This tutorial reproduces the **ideal (noiseless) simulation** of Fig. 4a from\n", + "[Bluvstein et al., Nature (2026)](https://www.nature.com/articles/s41586-025-09848-5):\n", + "subjecting quantum error-correction codes of different dimensionality to a\n", + "**global $R_Z(\\varphi)$ rotation**, then measuring the logical $X$ projection\n", + "and X-stabilizer revival.\n", + "\n", + "We use **Tsim** to prepare $|+_L\\rangle$ (hypercube encoder for $[[7,1,3]]$, Stim\n", + "graph-state synthesis for $[[15,1,3]]$). The rotation sweep uses exact state\n", + "vectors, which is fast for $n \\le 15$ and matches the paper's noiseless logical\n", + "response. It is **not** the experimental lookup-table decoder or postselection\n", + "pipeline." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75a69061", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from tsim import Circuit\n", + "\n", + "from utils.global_rotation import (\n", + " CODE_7_PLUS_L,\n", + " GATE_LABELS_DEG,\n", + " build_systems,\n", + " compute_curve,\n", + " reed_muller_prep_circuit,\n", + " verify_encoding_flows,\n", + " verify_steane_plus_l_flows,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "3a6179ad", + "metadata": {}, + "source": [ + "## Physical setup\n", + "\n", + "Fig. 4a studies\n", + "\n", + "$$|+_L\\rangle \\;\\xrightarrow{\\; R_Z(\\varphi)^{\\otimes n}\\;}\\; \\text{measure } X_L \\text{ and stabilizers}$$\n", + "\n", + "for four configurations:\n", + "\n", + "* **2D colour code** — $[[7,1,3]]$ Steane code\n", + "* **3D colour code** — $[[15,1,3]]$ quantum Reed–Muller code\n", + "* **Unentangled** — $|+\\rangle^{\\otimes 15}$ analysed with the same 3D observables\n", + "* **3D negative stabilizer** — Reed–Muller state with flipped Z-check signs\n", + "\n", + "Encoded 2D/3D codes show **X-stabilizer revivals** at Clifford angles ($0^\\circ,\n", + "\\pm 90^\\circ, \\ldots$). The 3D Reed–Muller code has an additional revival at\n", + "$\\pm 45^\\circ$ (transversal $T$), where $\\langle X_L\\rangle = 1/\\sqrt{2}$.\n", + "Unentangled qubits analysed with the same 15-body observables revive only near\n", + "$\\pm 180^\\circ$.\n", + "\n", + "This notebook is an **ideal noiseless** tutorial calculation. The experimental\n", + "Fig. 4a curves additionally use lookup-table decoding, postselection, and\n", + "purity normalization (Methods)." + ] + }, + { + "cell_type": "markdown", + "id": "c7c5d1e0", + "metadata": {}, + "source": [ + "## Hypercube encoding of $|+_L\\rangle$ for $[[7,1,3]]$\n", + "\n", + "The QuEra codes belong to the Reed–Muller family (Extended Data Fig. 10a). For\n", + "the 2D Steane code the **fourth CZ layer is turned off**. Qubit 6 starts in\n", + "$|+\\rangle$; qubits 0–5 in $|0\\rangle$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60e04d88", + "metadata": {}, + "outputs": [], + "source": [ + "encoding = Circuit(CODE_7_PLUS_L)\n", + "encoding.diagram(\"timeline-svg\", height=320)" + ] + }, + { + "cell_type": "markdown", + "id": "85462598", + "metadata": {}, + "source": [ + "Verify stabilizers and logical operators with Stim's `flow_generators()`, as\n", + "suggested in [issue #119](https://github.com/QuEraComputing/tsim/issues/119):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f088d17", + "metadata": {}, + "outputs": [], + "source": [ + "encoding_flows = verify_encoding_flows(CODE_7_PLUS_L)\n", + "prep_flows = verify_steane_plus_l_flows()\n", + "rm_flows = reed_muller_prep_circuit().flow_generators()\n", + "\n", + "print(f\"{len(encoding_flows)} flows on the Clifford encoding circuit\")\n", + "print(f\"{len(prep_flows)} flows on the full $[[7,1,3]]$ $|+_L\\\\rangle$ preparation\")\n", + "print(f\"{len(rm_flows)} flows on the $[[15,1,3]]$ preparation\")\n", + "for flow in encoding_flows[:4]:\n", + " print(flow)" + ] + }, + { + "cell_type": "markdown", + "id": "c5fed406", + "metadata": {}, + "source": [ + "## Global rotation sweep\n", + "\n", + "For each $\\varphi/\\pi \\in [-1, 1]$, we apply global $R_Z(\\varphi)$ and evaluate\n", + "\n", + "1. $\\langle X_L\\rangle$ (logical $X$ projection), and\n", + "2. mean $|\\langle S_j^X\\rangle|$ over X-type stabilizers (normalized to unit peak).\n", + "\n", + "Global $R_Z$ disturbs X-checks while leaving Z-checks invariant; the bottom panel\n", + "tracks X-stabilizer revivals, matching the Fig. 4a convention." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f351602", + "metadata": {}, + "outputs": [], + "source": [ + "N_ANGLES = 37\n", + "angles_deg = np.linspace(-180, 180, N_ANGLES)\n", + "\n", + "systems = build_systems()\n", + "curves = {system.label: compute_curve(system, angles_deg) for system in systems}" + ] + }, + { + "cell_type": "markdown", + "id": "214da9bd", + "metadata": {}, + "source": [ + "## Fig. 4a (ideal simulation)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e923c263", + "metadata": {}, + "outputs": [], + "source": [ + "plot_order = [\n", + " \"2D colour ($[[7,1,3]]$)\",\n", + " \"3D colour ($[[15,1,3]]$)\",\n", + " \"Unentangled\",\n", + " \"3D negative stabilizer\",\n", + "]\n", + "colors = [\"C0\", \"C3\", \"C2\", \"C1\"]\n", + "\n", + "fig, axes = plt.subplots(2, 1, figsize=(8, 5.5), sharex=True, constrained_layout=True)\n", + "\n", + "for label, color in zip(plot_order, colors):\n", + " curve = curves[label]\n", + " axes[0].plot(curve.angles_deg, curve.logical_x, \"-\", color=color, label=label)\n", + " if label != \"Unentangled\":\n", + " axes[1].plot(curve.angles_deg, curve.stabilizer_abs, \"-\", color=color, label=label)\n", + "\n", + "axes[0].set_ylabel(r\"Logical $X$ projection\")\n", + "axes[0].set_ylim(-1.05, 1.05)\n", + "axes[0].axhline(0, color=\"0.85\", lw=0.8)\n", + "axes[0].legend(loc=\"lower center\", ncol=2, fontsize=8, frameon=False)\n", + "\n", + "axes[1].set_ylabel(\"Normalized abs. stabilizer\")\n", + "axes[1].set_xlabel(r\"Global $\\varphi\\;(^{\\circ})$\")\n", + "axes[1].set_ylim(-0.05, 1.05)\n", + "axes[1].legend(loc=\"lower center\", ncol=2, fontsize=8, frameon=False)\n", + "\n", + "for ax in axes:\n", + " ax.set_xlim(-185, 185)\n", + " ax.set_xticks(list(GATE_LABELS_DEG.keys()))\n", + " ax.set_xticklabels(list(GATE_LABELS_DEG.values()))\n", + "\n", + "fig.suptitle(\"Fig. 4a — global $R_Z$ rotation (noiseless simulation)\", y=1.02)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b2938e84", + "metadata": {}, + "source": [ + "## What to notice\n", + "\n", + "* **2D colour** — X-stabilizers revive at $0^\\circ$ and $\\pm 90^\\circ$; logical\n", + " $X$ is $+1$ at $0^\\circ$ and $\\approx 0$ at $\\pm 90^\\circ$.\n", + "* **3D colour** — same Clifford stabilizer revivals **plus** a logical-$X$ value\n", + " $1/\\sqrt{2}$ at $\\pm 45^\\circ$ (transversal $T$).\n", + "* **Unentangled** — $\\langle X_L\\rangle = \\cos^{15}\\varphi$; revivals only near\n", + " $\\pm 180^\\circ$, not at Clifford angles.\n", + "* **Negative stabilizer** — wrong Z-check signs remove the 45° T feature while\n", + " Clifford stabilizer revivals remain.\n", + "\n", + "The $[[15,1,3]]$ state is prepared from verified stabilizer generators (Stim\n", + "graph-state synthesis). The paper's hypercube encoder (Extended Data Fig. 10a) is\n", + "shown and verified for $[[7,1,3]]$ above; the full 15-qubit hypercube circuit is\n", + "not in the supplementary Stim files." + ] + }, + { + "cell_type": "markdown", + "id": "494eba68", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "In this notebook we:\n", + "\n", + "* built and verified the QuEra hypercube encoder for $|+_L\\rangle$ on $[[7,1,3]]$,\n", + "* prepared $[[15,1,3]]$ $|+_L\\rangle$ from its stabilizer generators,\n", + "* swept a global $R_Z$ rotation and computed logical-$X$ and X-stabilizer observables, and\n", + "* reproduced the qualitative structure of Fig. 4a for all four curves." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.13.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/demos/utils/global_rotation.py b/docs/demos/utils/global_rotation.py new file mode 100644 index 00000000..ad5ddd49 --- /dev/null +++ b/docs/demos/utils/global_rotation.py @@ -0,0 +1,420 @@ +"""Fig. 4a global-rotation demo: ideal state-vector sweep for Reed–Muller codes. + +Tsim prepares $|+_L\\rangle$; global $R_Z(\\varphi)^{\\otimes n}$ is applied as a +diagonal phase on the state vector (exact, fast for $n \\le 15$). This matches the +paper's noiseless logical response; it is not the experimental decoder pipeline. +""" + +from __future__ import annotations + +import itertools +from dataclasses import dataclass +from functools import cache + +import numpy as np +import stim +import tsim + +# Hypercube encoder for $|+_L\\rangle$ on $[[7,1,3]]$ (Extended Data Fig. 10a; 4th CZ off). +CODE_7_PLUS_L = """ +RX 6 +R 0 1 2 3 4 5 +SQRT_Y_DAG 0 1 2 3 4 5 +CZ 1 2 3 4 5 6 +SQRT_Y 6 +CZ 0 3 2 5 4 6 +SQRT_Y 2 3 4 5 6 +CZ 0 1 2 3 4 5 +SQRT_Y 1 2 4 +""" + +STEANE_PREP = f"RX 6\n{CODE_7_PLUS_L.strip()}\n" + +# Same 4-weight supports as ``encoding_demo.ipynb`` (X and Z checks coincide). +STEANE_STABILIZER_SUPPORTS = [[0, 1, 2, 3], [1, 2, 4, 5], [2, 3, 4, 6]] +STEANE_LOGICAL_X_SUPPORT = [0, 1, 5] + +GATE_LABELS_DEG = { + -180: r"$Z$", + -135: r"$(TS)^\dagger$", + -90: r"$S^\dagger$", + -45: r"$T^\dagger$", + 0: r"$I$", + 45: r"$T$", + 90: r"$S$", + 135: r"$TS$", + 180: r"$Z$", +} + +PAULI_MUL: dict[tuple[str, str], tuple[complex, str]] = { + ("I", "I"): (1, "I"), + ("I", "X"): (1, "X"), + ("I", "Y"): (1, "Y"), + ("I", "Z"): (1, "Z"), + ("X", "I"): (1, "X"), + ("X", "X"): (1, "I"), + ("X", "Y"): (1j, "Z"), + ("X", "Z"): (-1j, "Y"), + ("Y", "I"): (1, "Y"), + ("Y", "X"): (-1j, "Z"), + ("Y", "Y"): (1, "I"), + ("Y", "Z"): (1j, "X"), + ("Z", "I"): (1, "Z"), + ("Z", "X"): (1j, "Y"), + ("Z", "Y"): (-1j, "X"), + ("Z", "Z"): (1, "I"), +} + + +@dataclass(frozen=True) +class Observable: + coefficient: complex + pauli: str + + +@dataclass(frozen=True) +class RotationSystem: + label: str + n: int + base_state: np.ndarray + x_stabilizers: list[str] + logical_x: Observable + logical_y: Observable + + +@dataclass(frozen=True) +class RotationCurve: + angles_deg: np.ndarray + logical_x: np.ndarray + stabilizer_abs: np.ndarray + + +def pauli_string(n: int, axis: str, support: list[int] | np.ndarray) -> str: + text = ["I"] * n + for qubit in support: + text[int(qubit)] = axis + return "".join(text) + + +def stim_pauli_text(pauli: str) -> str: + return pauli.replace("I", "_") + + +def multiply_paulis(left: str, right: str) -> tuple[complex, str]: + phase = 1 + out: list[str] = [] + for a, b in zip(left, right, strict=True): + local_phase, local_pauli = PAULI_MUL[(a, b)] + phase *= local_phase + out.append(local_pauli) + return phase, "".join(out) + + +def generated_group(generators: list[str]) -> list[tuple[complex, str]]: + identity = "I" * len(generators[0]) + group: list[tuple[complex, str]] = [] + for bits in itertools.product([0, 1], repeat=len(generators)): + phase = 1 + pauli = identity + for bit, generator in zip(bits, generators, strict=True): + if bit: + local_phase, pauli = multiply_paulis(pauli, generator) + phase *= local_phase + group.append((phase, pauli)) + return group + + +@cache +def pauli_action(pauli: str) -> tuple[np.ndarray, np.ndarray]: + n = len(pauli) + indices = np.arange(1 << n, dtype=np.uint64) + flip_mask = 0 + z_mask = 0 + y_count = 0 + for qubit, axis in enumerate(pauli): + bit_position = n - 1 - qubit + if axis in "XY": + flip_mask |= 1 << bit_position + if axis in "YZ": + z_mask |= 1 << bit_position + if axis == "Y": + y_count += 1 + parity = np.array( + [((int(index) & z_mask).bit_count() & 1) for index in indices], + dtype=np.int8, + ) + phase = ((1j) ** y_count) * (1 - 2 * parity) + return indices ^ np.uint64(flip_mask), phase + + +@cache +def hamming_weights(n: int) -> np.ndarray: + return np.array([index.bit_count() for index in range(1 << n)], dtype=np.int16) + + +def pauli_expectation_complex(state: np.ndarray, pauli: str) -> complex: + partner, phase = pauli_action(pauli) + return np.sum(np.conjugate(state[partner]) * phase * state) + + +def pauli_expectation(state: np.ndarray, pauli: str) -> float: + return float(np.real_if_close(pauli_expectation_complex(state, pauli)).real) + + +def observable_expectation(state: np.ndarray, observable: Observable) -> float: + value = observable.coefficient * pauli_expectation_complex(state, observable.pauli) + return float(np.real_if_close(value).real) + + +def apply_pauli_to_state(state: np.ndarray, pauli: str) -> np.ndarray: + partner, phase = pauli_action(pauli) + out = np.empty_like(state) + out[partner] = phase * state + return out + + +def logical_y_from_xz(logical_x: str, logical_z: str) -> Observable: + phase, pauli = multiply_paulis(logical_x, logical_z) + coefficient = 1j * phase + return Observable(float(np.real(coefficient)), pauli) + + +def normalized_state(circuit: tsim.Circuit) -> np.ndarray: + state = np.asarray(circuit.to_matrix()).reshape(-1) + return state / np.sqrt(np.vdot(state, state).real) + + +def rotate_state_from_base( + base_state: np.ndarray, n: int, angle_degrees: float +) -> np.ndarray: + theta = np.deg2rad(angle_degrees) + phases = np.exp(-0.5j * theta * (n - 2 * hamming_weights(n))) + return base_state * phases + + +def mean_abs_expectation(state: np.ndarray, observables: list[str]) -> float: + return float( + np.mean([abs(pauli_expectation(state, obs)) for obs in observables]) + ) + + +def normalize_signed(values: np.ndarray) -> np.ndarray: + scale = np.max(np.abs(values)) + if scale == 0: + return values + return values / scale + + +def verify_encoding_flows(encoding: str) -> list[stim.Flow]: + """Pauli flows for the Clifford hypercube encoding subcircuit.""" + return list(stim.Circuit(encoding).flow_generators()) + + +def verify_steane_plus_l_flows() -> list[stim.Flow]: + """Pauli flows for $|+_L\\rangle$ preparation on $[[7,1,3]]$.""" + return list(stim.Circuit(STEANE_PREP).flow_generators()) + + +def gf2_rank(rows: list[np.ndarray] | np.ndarray) -> int: + rows_array = np.asarray(rows, dtype=np.uint8) + if rows_array.ndim == 1: + rows_array = rows_array.reshape(1, -1) + n = rows_array.shape[1] + rows_int = [ + int(sum(int(bit) << i for i, bit in enumerate(row))) + for row in rows_array + if np.any(row) + ] + rank = 0 + for column in range(n): + pivot = next( + (i for i in range(rank, len(rows_int)) if (rows_int[i] >> column) & 1), + None, + ) + if pivot is None: + continue + rows_int[rank], rows_int[pivot] = rows_int[pivot], rows_int[rank] + for i in range(len(rows_int)): + if i != rank and ((rows_int[i] >> column) & 1): + rows_int[i] ^= rows_int[rank] + rank += 1 + return rank + + +def gf2_nullspace_basis(matrix: np.ndarray) -> list[np.ndarray]: + matrix = np.asarray(matrix, dtype=np.uint8) + n = matrix.shape[1] + rows = [int(sum(int(bit) << i for i, bit in enumerate(row))) for row in matrix] + pivots: list[int] = [] + row_index = 0 + for column in range(n): + pivot = next( + (i for i in range(row_index, len(rows)) if (rows[i] >> column) & 1), + None, + ) + if pivot is None: + continue + rows[row_index], rows[pivot] = rows[pivot], rows[row_index] + for i in range(len(rows)): + if i != row_index and ((rows[i] >> column) & 1): + rows[i] ^= rows[row_index] + pivots.append(column) + row_index += 1 + + basis: list[np.ndarray] = [] + for free_column in (c for c in range(n) if c not in pivots): + vector = 1 << free_column + for row_number, pivot_column in enumerate(pivots): + if (rows[row_number] >> free_column) & 1: + vector |= 1 << pivot_column + basis.append(np.array([(vector >> i) & 1 for i in range(n)], dtype=np.uint8)) + return basis + + +def build_reed_muller_15() -> tuple[ + list[str], + list[str], + Observable, + Observable, + np.ndarray, +]: + """Return X/Z stabilizers, logical X/Y, and $|+_L\\rangle$ state.""" + rm_n = 15 + rm_columns = np.array( + [[int(bit) for bit in f"{value:04b}"] for value in range(1, 16)], + dtype=np.uint8, + ).T + rm_logical_x_row = np.ones(rm_n, dtype=np.uint8) + rm_x_rows = [row.copy() for row in rm_columns] + rm_z_rows = gf2_nullspace_basis(np.vstack([rm_logical_x_row, *rm_x_rows])) + + rm_logical_z_row: np.ndarray | None = None + for support in itertools.combinations(range(rm_n), 3): + candidate = np.zeros(rm_n, dtype=np.uint8) + candidate[list(support)] = 1 + if all(np.dot(candidate, row) % 2 == 0 for row in rm_x_rows) and ( + np.dot(candidate, rm_logical_x_row) % 2 == 1 + ): + rm_logical_z_row = candidate + break + assert rm_logical_z_row is not None + + x_stabilizers = [ + pauli_string(rm_n, "X", np.flatnonzero(row)) for row in rm_x_rows + ] + z_stabilizers = [ + pauli_string(rm_n, "Z", np.flatnonzero(row)) for row in rm_z_rows + ] + logical_x = Observable( + 1, pauli_string(rm_n, "X", np.flatnonzero(rm_logical_x_row)) + ) + logical_z = pauli_string(rm_n, "Z", np.flatnonzero(rm_logical_z_row)) + + stabilizer_strings = x_stabilizers + z_stabilizers + [logical_x.pauli] + stim_stabilizers = [ + stim.PauliString("+" + stim_pauli_text(p)) for p in stabilizer_strings + ] + prep_stim = stim.Tableau.from_stabilizers( + stim_stabilizers, + allow_underconstrained=False, + allow_redundant=False, + ).to_circuit(method="graph_state") + prep_circuit = tsim.Circuit.from_stim_program(prep_stim) + base_state = normalized_state(prep_circuit) + + logical_y = logical_y_from_xz(logical_x.pauli, logical_z) + test_state = rotate_state_from_base(base_state, rm_n, 1.0) + if observable_expectation(test_state, logical_y) < 0: + logical_y = Observable(-logical_y.coefficient, logical_y.pauli) + + return x_stabilizers, z_stabilizers, logical_x, logical_y, base_state + + +def reed_muller_prep_circuit() -> stim.Circuit: + """Stim preparation circuit for $|+_L\\rangle$ on $[[15,1,3]]$.""" + rm_x, rm_z, logical_x, _, _ = build_reed_muller_15() + stabilizers = rm_x + rm_z + [logical_x.pauli] + return stim.Tableau.from_stabilizers( + [stim.PauliString("+" + stim_pauli_text(p)) for p in stabilizers], + allow_underconstrained=False, + allow_redundant=False, + ).to_circuit(method="graph_state") + + +def build_steane_7() -> tuple[list[str], Observable, Observable, np.ndarray]: + """Hypercube $|+_L\\rangle$ on $[[7,1,3]]$ with Steane stabilizer labelling.""" + n = 7 + x_stabilizers = [ + pauli_string(n, "X", support) for support in STEANE_STABILIZER_SUPPORTS + ] + logical_x = Observable(1, pauli_string(n, "X", STEANE_LOGICAL_X_SUPPORT)) + logical_z = pauli_string(n, "Z", STEANE_LOGICAL_X_SUPPORT) + logical_y = logical_y_from_xz(logical_x.pauli, logical_z) + + base_state = normalized_state(tsim.Circuit(STEANE_PREP)) + test_state = rotate_state_from_base(base_state, n, 1.0) + if observable_expectation(test_state, logical_y) < 0: + logical_y = Observable(-logical_y.coefficient, logical_y.pauli) + + return x_stabilizers, logical_x, logical_y, base_state + + +def build_systems() -> list[RotationSystem]: + steane_x, steane_lx, steane_ly, steane_state = build_steane_7() + rm_x, _, rm_lx, rm_ly, rm_state = build_reed_muller_15() + + product_state = normalized_state( + tsim.Circuit("RX " + " ".join(map(str, range(15)))) + ) + + rm_corner_flip = pauli_string(15, "X", [0, 1, 3, 7]) + rm_negative_state = apply_pauli_to_state(rm_state, rm_corner_flip) + + return [ + RotationSystem( + "2D colour ($[[7,1,3]]$)", + 7, + steane_state, + steane_x, + steane_lx, + steane_ly, + ), + RotationSystem( + "3D colour ($[[15,1,3]]$)", + 15, + rm_state, + rm_x, + rm_lx, + rm_ly, + ), + RotationSystem( + "Unentangled", + 15, + product_state, + rm_x, + rm_lx, + rm_ly, + ), + RotationSystem( + "3D negative stabilizer", + 15, + rm_negative_state, + rm_x, + rm_lx, + rm_ly, + ), + ] + + +def compute_curve(system: RotationSystem, angles_deg: np.ndarray) -> RotationCurve: + logical_x_vals: list[float] = [] + stab_vals: list[float] = [] + for angle in angles_deg: + state = rotate_state_from_base(system.base_state, system.n, float(angle)) + logical_x_vals.append(observable_expectation(state, system.logical_x)) + stab_vals.append(mean_abs_expectation(state, system.x_stabilizers)) + return RotationCurve( + angles_deg=np.asarray(angles_deg, dtype=float), + logical_x=normalize_signed(np.asarray(logical_x_vals)), + stabilizer_abs=normalize_signed(np.asarray(stab_vals)), + ) diff --git a/mkdocs.yml b/mkdocs.yml index acc5a219..8e012fbb 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -15,6 +15,7 @@ nav: - Contributing: contrib.md - Tutorials: - demos/encoding_demo.ipynb + - demos/global_rotations.ipynb - demos/magic_state_distillation.ipynb - demos/from_stim_to_tsim.ipynb - API Reference: reference/ diff --git a/test/unit/demos/test_global_rotation.py b/test/unit/demos/test_global_rotation.py new file mode 100644 index 00000000..83112217 --- /dev/null +++ b/test/unit/demos/test_global_rotation.py @@ -0,0 +1,45 @@ +"""Smoke tests for the Fig. 4a global-rotation demo utilities.""" + +import sys +from pathlib import Path + +import numpy as np + +_DEMOS = Path(__file__).resolve().parents[3] / "docs" / "demos" +sys.path.insert(0, str(_DEMOS)) + +from utils.global_rotation import ( # noqa: E402 + CODE_7_PLUS_L, + build_systems, + compute_curve, + verify_encoding_flows, + verify_steane_plus_l_flows, +) + +KEY_ANGLES = np.array([-180, -90, -45, 0, 45, 90, 180], dtype=float) + + +def test_steane_flow_counts(): + assert len(verify_encoding_flows(CODE_7_PLUS_L)) == 7 + assert len(verify_steane_plus_l_flows()) == 7 + + +def test_fig4a_qualitative_curves(): + systems = {s.label: s for s in build_systems()} + curves = {label: compute_curve(systems[label], KEY_ANGLES) for label in systems} + + xl_2d = curves["2D colour ($[[7,1,3]]$)"].logical_x + xl_3d = curves["3D colour ($[[15,1,3]]$)"].logical_x + xl_u = curves["Unentangled"].logical_x + stab_2d = curves["2D colour ($[[7,1,3]]$)"].stabilizer_abs + stab_3d = curves["3D colour ($[[15,1,3]]$)"].stabilizer_abs + + assert abs(xl_2d[3] - 1.0) < 0.02 + assert abs(xl_3d[3] - 1.0) < 0.02 + assert abs(xl_3d[4] - 0.707) < 0.05 + assert abs(xl_u[0] + 1.0) < 0.02 + assert abs(xl_u[3] - 1.0) < 0.02 + assert stab_2d[5] > 0.9 + assert stab_3d[5] > 0.9 + assert stab_3d[4] > 0.9 + assert xl_u[5] < 0.1 From b8943f590cb972fa3420acf82412f0782ef58131 Mon Sep 17 00:00:00 2001 From: rosspeili Date: Tue, 9 Jun 2026 16:53:59 +0300 Subject: [PATCH 2/2] Guard normalized_state against Clifford-only to_matrix flattening. --- docs/demos/utils/global_rotation.py | 45 +++++++++++++++++++++++++---- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/docs/demos/utils/global_rotation.py b/docs/demos/utils/global_rotation.py index ad5ddd49..6b75b95d 100644 --- a/docs/demos/utils/global_rotation.py +++ b/docs/demos/utils/global_rotation.py @@ -179,8 +179,41 @@ def logical_y_from_xz(logical_x: str, logical_z: str) -> Observable: return Observable(float(np.real(coefficient)), pauli) -def normalized_state(circuit: tsim.Circuit) -> np.ndarray: - state = np.asarray(circuit.to_matrix()).reshape(-1) +_STATE_INIT_OPS = frozenset( + {"R", "RX", "RY", "RZ", "MR", "MRX", "MRY", "MRZ", "MPP", "MXX", "MYY", "MZZ"} +) + + +def _has_state_initialization(stim_circuit: stim.Circuit) -> bool: + return any(instruction.name in _STATE_INIT_OPS for instruction in stim_circuit) + + +def normalized_state(circuit: tsim.Circuit, *, n_qubits: int | None = None) -> np.ndarray: + """Return a normalized state vector from a Tsim preparation circuit. + + ``Circuit.to_matrix()`` returns a column state vector only when the Stim + program includes initialization (``R`` / ``RX`` / …). A Clifford-only circuit + would materialize the full unitary and ``reshape(-1)`` would yield length + ``2**(2n)``, which is both expensive and incompatible with our phase sweep. + """ + stim_circuit = circuit.stim_circuit + if not _has_state_initialization(stim_circuit): + raise ValueError( + "Preparation circuit has no R/RX-style initialization; " + "cannot extract a state vector via to_matrix()." + ) + + matrix = np.asarray(circuit.to_matrix()) + if matrix.ndim == 2 and matrix.shape[0] == matrix.shape[1] and matrix.shape[1] != 1: + raise ValueError( + f"Expected a state vector from to_matrix(), got unitary shape " + f"{matrix.shape}." + ) + state = matrix.reshape(-1) + if n_qubits is not None and state.size != 1 << n_qubits: + raise ValueError( + f"Expected state vector of length {1 << n_qubits}, got {state.size}." + ) return state / np.sqrt(np.vdot(state, state).real) @@ -319,8 +352,9 @@ def build_reed_muller_15() -> tuple[ allow_underconstrained=False, allow_redundant=False, ).to_circuit(method="graph_state") + # graph_state synthesis prepends RX on all qubits before the Clifford layer. prep_circuit = tsim.Circuit.from_stim_program(prep_stim) - base_state = normalized_state(prep_circuit) + base_state = normalized_state(prep_circuit, n_qubits=rm_n) logical_y = logical_y_from_xz(logical_x.pauli, logical_z) test_state = rotate_state_from_base(base_state, rm_n, 1.0) @@ -351,7 +385,7 @@ def build_steane_7() -> tuple[list[str], Observable, Observable, np.ndarray]: logical_z = pauli_string(n, "Z", STEANE_LOGICAL_X_SUPPORT) logical_y = logical_y_from_xz(logical_x.pauli, logical_z) - base_state = normalized_state(tsim.Circuit(STEANE_PREP)) + base_state = normalized_state(tsim.Circuit(STEANE_PREP), n_qubits=n) test_state = rotate_state_from_base(base_state, n, 1.0) if observable_expectation(test_state, logical_y) < 0: logical_y = Observable(-logical_y.coefficient, logical_y.pauli) @@ -364,7 +398,8 @@ def build_systems() -> list[RotationSystem]: rm_x, _, rm_lx, rm_ly, rm_state = build_reed_muller_15() product_state = normalized_state( - tsim.Circuit("RX " + " ".join(map(str, range(15)))) + tsim.Circuit("RX " + " ".join(map(str, range(15)))), + n_qubits=15, ) rm_corner_flip = pauli_string(15, "X", [0, 1, 3, 7])