diff --git a/CHANGELOG.md b/CHANGELOG.md index 44cabffa..831f622a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,12 +10,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - **QEC Stim MPP Import**: Added utilities for building `StabilizerCode` inputs from unsigned Stim `MPP` layers, including sparse Stim qubit id mapping, coordinate import, multi-layer selection, detector/logical-observable import, and the optional `graphqomb[stim]` extra. Signed products using inverted Pauli targets are rejected because stabilizer signs are not retained. +- **Input Initialization Bases**: Added per-input positive Pauli eigenstate initialization via `GraphState.register_input(..., init_axis=...)` and `assign_input_initialization_axis()`, with `X+`, `Y+`, and `Z+` propagated through `qompile()`, `Pattern`, `PatternSimulator`, and Stim export. ### Changed - **Development Tooling**: Use uv as the default dependency manager for local development, CI, documentation builds, and publishing workflows. - **Graph State API**: Replaced legacy `physical_*` graph methods/properties with standard graph-style `nodes`, `edges`, `add_node()`, `add_edge()`, `remove_node()`, `remove_edge()`, and count/query helpers. - **Pattern Simulator**: Materialize pending output Pauli-frame corrections when returning output statevectors or explicit output measurement results. +- **PTN Format**: Bumped exported `.ptn` files to format version 2 to record non-default input initialization bases with `.input_basis`; version 1 files remain readable and default inputs to `X+`. +- **Stim Compiler**: Input reset instructions now follow the stored initialization basis (`RX`, `RY`, or `R`) while non-input preparations remain `RX`. ### Removed diff --git a/graphqomb/graphstate.py b/graphqomb/graphstate.py index 37c5853a..4c29ec0d 100644 --- a/graphqomb/graphstate.py +++ b/graphqomb/graphstate.py @@ -26,7 +26,7 @@ import typing_extensions -from graphqomb.common import MeasBasis, Plane, PlannerMeasBasis +from graphqomb.common import Axis, MeasBasis, Plane, PlannerMeasBasis from graphqomb.euler import update_lc_basis, update_lc_lc if TYPE_CHECKING: @@ -49,6 +49,17 @@ def input_node_indices(self) -> dict[int, int]: qubit indices map of input nodes. """ + @property + def input_initialization_axes(self) -> dict[int, Axis]: + r"""Input initialization Pauli axes. + + Returns + ------- + `dict`\[`int`, `Axis`\] + map of input nodes to Pauli initialization axes. + """ + return dict.fromkeys(self.input_node_indices, Axis.X) + @property @abc.abstractmethod def output_node_indices(self) -> dict[int, int]: @@ -174,7 +185,7 @@ def number_of_edges(self) -> int: return len(self.edges) @abc.abstractmethod - def register_input(self, node: int, q_index: int) -> None: + def register_input(self, node: int, q_index: int, *, init_axis: Axis = Axis.X) -> None: """Mark the node as an input node. Parameters @@ -183,6 +194,8 @@ def register_input(self, node: int, q_index: int) -> None: node index q_index : `int` logical qubit index + init_axis : `Axis`, optional + Pauli axis for positive-eigenstate initialization, by default Axis.X """ @abc.abstractmethod @@ -244,6 +257,7 @@ class GraphState(BaseGraphState): """Minimal implementation of GraphState.""" __input_node_indices: dict[int, int] + __input_initialization_axes: dict[int, Axis] __output_node_indices: dict[int, int] __nodes: set[int] __neighbors: dict[int, set[int]] @@ -257,6 +271,7 @@ class GraphState(BaseGraphState): def __init__(self) -> None: self.__input_node_indices = {} + self.__input_initialization_axes = {} self.__output_node_indices = {} self.__nodes = set() self.__neighbors = {} @@ -278,6 +293,18 @@ def input_node_indices(self) -> dict[int, int]: """ return self.__input_node_indices.copy() + @property + @typing_extensions.override + def input_initialization_axes(self) -> dict[int, Axis]: + r"""Input initialization Pauli axes. + + Returns + ------- + `dict`\[`int`, `Axis`\] + map of input nodes to Pauli initialization axes. + """ + return {node: self.__input_initialization_axes.get(node, Axis.X) for node in self.__input_node_indices} + @property @typing_extensions.override def output_node_indices(self) -> dict[int, int]: @@ -483,6 +510,7 @@ def remove_node(self, node: int) -> None: if node in self.output_node_indices: del self.__output_node_indices[node] + self.__input_initialization_axes.pop(node, None) self.__meas_bases.pop(node, None) self.__local_cliffords.pop(node, None) self.__coordinates.pop(node, None) @@ -514,7 +542,7 @@ def remove_edge(self, node1: int, node2: int) -> None: self.__neighbors[node2] -= {node1} @typing_extensions.override - def register_input(self, node: int, q_index: int) -> None: + def register_input(self, node: int, q_index: int, *, init_axis: Axis = Axis.X) -> None: """Mark the node as an input node. Parameters @@ -523,6 +551,8 @@ def register_input(self, node: int, q_index: int) -> None: node index q_index : `int` logical qubit index + init_axis : `Axis`, optional + Pauli axis for positive-eigenstate initialization, by default Axis.X Raises ------ @@ -536,7 +566,44 @@ def register_input(self, node: int, q_index: int) -> None: if q_index in self.input_node_indices.values(): msg = "The q_index already exists in input qubit indices" raise ValueError(msg) + self._check_input_initialization_axis(init_axis) self.__input_node_indices[node] = q_index + self.__input_initialization_axes[node] = init_axis + + def assign_input_initialization_axis(self, node: int, axis: Axis) -> None: + """Assign a Pauli initialization axis to an input node. + + Parameters + ---------- + node : `int` + input node index + axis : `Axis` + Pauli axis for positive-eigenstate initialization + + Raises + ------ + ValueError + If the node does not exist or is not registered as an input. + """ + self._ensure_node_exists(node) + if node not in self.__input_node_indices: + msg = f"Node is not registered as an input node: {node}" + raise ValueError(msg) + self._check_input_initialization_axis(axis) + self.__input_initialization_axes[node] = axis + + @staticmethod + def _check_input_initialization_axis(axis: Axis) -> None: + """Validate an input initialization axis. + + Raises + ------ + TypeError + If the axis is not an Axis value. + """ + if not isinstance(axis, Axis): + msg = "Input initialization axis must be one of Axis.X, Axis.Y, Axis.Z" + raise TypeError(msg) @typing_extensions.override def register_output(self, node: int, q_index: int) -> None: @@ -677,14 +744,18 @@ def _expand_input_local_cliffords(self) -> dict[int, LocalCliffordExpansion]: """ node_index_addition_map: dict[int, LocalCliffordExpansion] = {} new_input_indices: dict[int, int] = {} + new_input_initialization_axes: dict[int, Axis] = {} for input_node, q_index in self.input_node_indices.items(): + init_axis = self.input_initialization_axes[input_node] lc = self._pop_local_clifford(input_node) if lc is None: new_input_indices[input_node] = q_index + new_input_initialization_axes[input_node] = init_axis continue new_node_index0 = self.add_node() new_input_indices[new_node_index0] = q_index + new_input_initialization_axes[new_node_index0] = init_axis new_node_index1 = self.add_node() new_node_index2 = self.add_node() @@ -701,8 +772,13 @@ def _expand_input_local_cliffords(self) -> dict[int, LocalCliffordExpansion]: ) self.__input_node_indices = {} + self.__input_initialization_axes = {} for new_input_index, q_index in new_input_indices.items(): - self.register_input(new_input_index, q_index) + self.register_input( + new_input_index, + q_index, + init_axis=new_input_initialization_axes[new_input_index], + ) return node_index_addition_map @@ -754,6 +830,7 @@ def from_graph( # noqa: C901, PLR0912, PLR0913, PLR0917 outputs: Sequence[NodeT] | None = None, meas_bases: Mapping[NodeT, MeasBasis] | None = None, coordinates: Mapping[NodeT, tuple[float, ...]] | None = None, + input_initialization_axes: Mapping[NodeT, Axis] | None = None, ) -> tuple[GraphState, dict[NodeT, int]]: r"""Create a graph state from nodes and edges with arbitrary node types. @@ -778,6 +855,8 @@ def from_graph( # noqa: C901, PLR0912, PLR0913, PLR0917 Default is None (no bases assigned initially). coordinates : `collections.abc.Mapping`\[NodeT, `tuple`\[`float`, ...\]\] | `None`, optional Coordinates for nodes (2D or 3D). Default is None (no coordinates). + input_initialization_axes : `collections.abc.Mapping`\[NodeT, `Axis`\] | `None`, optional + Pauli initialization axes for input nodes. Default is None (all inputs use Axis.X). Returns ------- @@ -844,7 +923,10 @@ def from_graph( # noqa: C901, PLR0912, PLR0913, PLR0917 # Register inputs with sequential qubit indices if inputs is not None: for q_index, input_node in enumerate(inputs): - graph_state.register_input(node_map[input_node], q_index) + init_axis = ( + Axis.X if input_initialization_axes is None else input_initialization_axes.get(input_node, Axis.X) + ) + graph_state.register_input(node_map[input_node], q_index, init_axis=init_axis) # Register outputs with sequential qubit indices if outputs is not None: @@ -908,7 +990,11 @@ def from_base_graph_state( # Register inputs with same qubit indices for input_node, q_index in base.input_node_indices.items(): - graph_state.register_input(node_map[input_node], q_index) + graph_state.register_input( + node_map[input_node], + q_index, + init_axis=base.input_initialization_axes.get(input_node, Axis.X), + ) # Register outputs with same qubit indices for output_node, q_index in base.output_node_indices.items(): @@ -1018,11 +1104,19 @@ def compose( # noqa: C901, PLR0912 # Register input nodes with preserved qindices for input_node, q_index in graph1.input_node_indices.items(): - composed_graph.register_input(node_map1[input_node], q_index) + composed_graph.register_input( + node_map1[input_node], + q_index, + init_axis=graph1.input_initialization_axes.get(input_node, Axis.X), + ) for input_node, q_index in graph2.input_node_indices.items(): if q_index not in target_q_indices: - composed_graph.register_input(node_map2[input_node], q_index) + composed_graph.register_input( + node_map2[input_node], + q_index, + init_axis=graph2.input_initialization_axes.get(input_node, Axis.X), + ) # Register output nodes with preserved qindices for output_node, q_index in graph1.output_node_indices.items(): diff --git a/graphqomb/pattern.py b/graphqomb/pattern.py index ef36f5d8..0063a302 100644 --- a/graphqomb/pattern.py +++ b/graphqomb/pattern.py @@ -20,9 +20,32 @@ if TYPE_CHECKING: from collections.abc import Callable, Iterator + from graphqomb.common import Axis from graphqomb.pauli_frame import PauliFrame +def _empty_input_coordinates() -> dict[int, tuple[float, ...]]: + r"""Return an empty input-coordinate map. + + Returns + ------- + `dict`\[`int`, `tuple`\[`float`, ...\]\] + Empty input-coordinate map. + """ + return {} + + +def _empty_input_initialization_axes() -> dict[int, Axis]: + r"""Return an empty input-initialization axis map. + + Returns + ------- + `dict`\[`int`, `Axis`\] + Empty input-initialization axis map. + """ + return {} + + @dataclasses.dataclass(frozen=True) class Pattern(Sequence[Command]): r"""Pattern class. @@ -39,13 +62,16 @@ class Pattern(Sequence[Command]): Pauli frame of the pattern to track the Pauli state of each node input_coordinates : `dict`\[`int`, `tuple`\[`float`, ...\]\] Coordinates for input nodes (2D or 3D) + input_initialization_axes : `dict`\[`int`, `Axis`\] + Pauli initialization axes for input nodes. Missing inputs default to Axis.X. """ input_node_indices: dict[int, int] output_node_indices: dict[int, int] commands: tuple[Command, ...] pauli_frame: PauliFrame - input_coordinates: dict[int, tuple[float, ...]] = dataclasses.field(default_factory=dict) + input_coordinates: dict[int, tuple[float, ...]] = dataclasses.field(default_factory=_empty_input_coordinates) + input_initialization_axes: dict[int, Axis] = dataclasses.field(default_factory=_empty_input_initialization_axes) def __len__(self) -> int: return len(self.commands) diff --git a/graphqomb/ptn_format.py b/graphqomb/ptn_format.py index 8765e667..76b38ed7 100644 --- a/graphqomb/ptn_format.py +++ b/graphqomb/ptn_format.py @@ -38,7 +38,8 @@ if TYPE_CHECKING: from collections.abc import Sequence -PTN_VERSION = 1 +PTN_VERSION = 2 +SUPPORTED_PTN_VERSIONS = frozenset({1, PTN_VERSION}) # Angle formatting/parsing lookup tables _ANGLE_TO_STR: dict[float, str] = { @@ -178,6 +179,13 @@ def _write_header(out: StringIO, pattern: Pattern) -> None: f"{node}:{qidx}" for node, qidx in sorted(pattern.input_node_indices.items(), key=operator.itemgetter(1)) ] out.write(f".input {' '.join(input_parts)}\n") + input_basis_parts = [ + f"{node}:{axis.name}" + for node, qidx in sorted(pattern.input_node_indices.items(), key=operator.itemgetter(1)) + if (axis := pattern.input_initialization_axes.get(node, Axis.X)) is not Axis.X + ] + if input_basis_parts: + out.write(f".input_basis {' '.join(input_basis_parts)}\n") if pattern.output_node_indices: output_parts = [ @@ -381,6 +389,38 @@ def _parse_node_qubit_pairs(parts: Sequence[str]) -> dict[int, int]: return result +def _parse_node_axis_pairs(parts: Sequence[str]) -> dict[int, Axis]: + r"""Parse node:axis pairs from string parts. + + Returns + ------- + `dict`\[`int`, `Axis`\] + Mapping from node to Pauli axis. + + Raises + ------ + ValueError + If any pair is malformed, duplicated, or uses an invalid axis. + """ + result: dict[int, Axis] = {} + for part in parts: + pair = part.split(":") + if len(pair) != 2: # noqa: PLR2004 + msg = f"Invalid node:axis pair: {part!r}" + raise ValueError(msg) + node_str, axis_str = pair + node = _parse_int(node_str, "node") + if node in result: + msg = f"Duplicate input basis node: {node}" + raise ValueError(msg) + try: + result[node] = Axis[axis_str] + except KeyError as exc: + msg = f"Invalid input basis axis: {axis_str!r}" + raise ValueError(msg) from exc + return result + + def _parse_node_set(parts: Sequence[str], label: str) -> set[int]: r"""Parse a non-empty set of node ids. @@ -454,6 +494,17 @@ def _empty_coordinates() -> dict[int, tuple[float, ...]]: return {} +def _empty_input_initialization_axes() -> dict[int, Axis]: + r"""Return an empty input-initialization axis map. + + Returns + ------- + `dict`\[`int`, `Axis`\] + Empty input-initialization axis map. + """ + return {} + + def _empty_commands() -> list[Command]: r"""Return an empty command list. @@ -499,6 +550,8 @@ class _PatternData: Mapping from node to qubit index for output nodes. input_coordinates : `dict`[`int`, `tuple`[`float`, ...]] Coordinates for input nodes. + input_initialization_axes : `dict`[`int`, `Axis`] + Pauli initialization axes for input nodes. commands : `list`[`Command`] List of quantum commands. xflow : `dict`[`int`, `set`[`int`]] @@ -512,6 +565,7 @@ class _PatternData: input_node_indices: dict[int, int] = field(default_factory=_empty_node_index_map) output_node_indices: dict[int, int] = field(default_factory=_empty_node_index_map) input_coordinates: dict[int, tuple[float, ...]] = field(default_factory=_empty_coordinates) + input_initialization_axes: dict[int, Axis] = field(default_factory=_empty_input_initialization_axes) commands: list[Command] = field(default_factory=_empty_commands) xflow: dict[int, set[int]] = field(default_factory=_empty_node_set_map) zflow: dict[int, set[int]] = field(default_factory=_empty_node_set_map) @@ -529,6 +583,7 @@ class _LoadedGraphState(BaseGraphState): _edges: set[tuple[int, int]] _meas_bases: dict[int, MeasBasis] _coordinates: dict[int, tuple[float, ...]] + _input_initialization_axes: dict[int, Axis] _neighbors: dict[int, set[int]] = field(init=False, repr=False) def __post_init__(self) -> None: @@ -538,6 +593,7 @@ def __post_init__(self) -> None: self._edges = {(node1, node2) if node1 < node2 else (node2, node1) for node1, node2 in self._edges} self._meas_bases = dict(self._meas_bases) self._coordinates = dict(self._coordinates) + self._input_initialization_axes = dict(self._input_initialization_axes) self._neighbors: dict[int, set[int]] = {node: set() for node in self._nodes} for node1, node2 in self._edges: self._neighbors.setdefault(node1, set()).add(node2) @@ -547,6 +603,10 @@ def __post_init__(self) -> None: def input_node_indices(self) -> dict[int, int]: return self._input_node_indices.copy() + @property + def input_initialization_axes(self) -> dict[int, Axis]: + return {node: self._input_initialization_axes.get(node, Axis.X) for node in self._input_node_indices} + @property def output_node_indices(self) -> dict[int, int]: return self._output_node_indices.copy() @@ -575,7 +635,7 @@ def add_edge(self, node1: int, node2: int) -> None: msg = "Loaded .ptn graph states are read-only" raise NotImplementedError(msg) - def register_input(self, node: int, q_index: int) -> None: + def register_input(self, node: int, q_index: int, *, init_axis: Axis = Axis.X) -> None: msg = "Loaded .ptn graph states are read-only" raise NotImplementedError(msg) @@ -615,6 +675,27 @@ def _command_nodes(cmd: Command) -> set[int]: return set() +def _input_initialization_axes_from_data(data: _PatternData) -> dict[int, Axis]: + r"""Validate and normalize parsed input initialization axes. + + Returns + ------- + `dict`\[`int`, `Axis`\] + Normalized input initialization axes. + + Raises + ------ + ValueError + If an input basis is specified for a non-input node. + """ + non_input_basis_nodes = set(data.input_initialization_axes) - set(data.input_node_indices) + if non_input_basis_nodes: + msg = f"Input basis specified for non-input node(s): {sorted(non_input_basis_nodes)}" + raise ValueError(msg) + + return {node: data.input_initialization_axes.get(node, Axis.X) for node in data.input_node_indices} + + def _build_pattern(data: _PatternData) -> Pattern: """Build a Pattern from parsed .ptn data. @@ -628,6 +709,8 @@ def _build_pattern(data: _PatternData) -> Pattern: ValueError If parsed commands contain invalid graph structure. """ + input_initialization_axes = _input_initialization_axes_from_data(data) + nodes: set[int] = set(data.input_node_indices) | set(data.output_node_indices) | set(data.input_coordinates) edges: set[tuple[int, int]] = set() meas_bases: dict[int, MeasBasis] = {} @@ -665,6 +748,7 @@ def _build_pattern(data: _PatternData) -> Pattern: _edges=edges, _meas_bases=meas_bases, _coordinates=coordinates, + _input_initialization_axes=input_initialization_axes, ) pauli_frame = PauliFrame( graphstate, @@ -679,6 +763,7 @@ def _build_pattern(data: _PatternData) -> Pattern: commands=tuple(data.commands), pauli_frame=pauli_frame, input_coordinates=dict(data.input_coordinates), + input_initialization_axes=input_initialization_axes, ) @@ -688,7 +773,7 @@ class _Parser: def __init__(self) -> None: self.result = _PatternData() self.current_timeslice = -1 - self.version_found = False + self.version: int | None = None def parse(self, s: str) -> Pattern: r"""Parse the input string and return Pattern. @@ -711,9 +796,12 @@ def parse(self, s: str) -> Pattern: for line_num, raw_line in enumerate(s.splitlines(), 1): self._parse_line(line_num, raw_line) - if not self.version_found: + if self.version is None: msg = "Missing .version directive" raise ValueError(msg) + if self.version == 1 and self.result.input_initialization_axes: + msg = ".input_basis requires .ptn version 2 or later" + raise ValueError(msg) return _build_pattern(self.result) @@ -759,6 +847,8 @@ def _parse_directive(self, line: str) -> None: self._handle_version(content) elif directive == ".input": self.result.input_node_indices = _parse_node_qubit_pairs(content.split()) + elif directive == ".input_basis": + self.result.input_initialization_axes = _parse_node_axis_pairs(content.split()) elif directive == ".output": self.result.output_node_indices = _parse_node_qubit_pairs(content.split()) elif directive == ".coord": @@ -787,10 +877,11 @@ def _handle_version(self, content: str) -> None: If the version is unsupported. """ version = _parse_int(content, "version") - if version != PTN_VERSION: - msg = f"Unsupported .ptn version: {version} (expected {PTN_VERSION})" + if version not in SUPPORTED_PTN_VERSIONS: + supported = ", ".join(str(supported_version) for supported_version in sorted(SUPPORTED_PTN_VERSIONS)) + msg = f"Unsupported .ptn version: {version} (supported: {supported})" raise ValueError(msg) - self.version_found = True + self.version = version def _handle_coord(self, content: str) -> None: """Handle .coord directive. diff --git a/graphqomb/qompiler.py b/graphqomb/qompiler.py index d49dec39..14858a92 100644 --- a/graphqomb/qompiler.py +++ b/graphqomb/qompiler.py @@ -149,4 +149,5 @@ def _qompile( commands=tuple(commands), pauli_frame=pauli_frame, input_coordinates=input_coords, + input_initialization_axes=graph.input_initialization_axes, ) diff --git a/graphqomb/simulator.py b/graphqomb/simulator.py index 559bc2ec..7ef0da20 100644 --- a/graphqomb/simulator.py +++ b/graphqomb/simulator.py @@ -16,13 +16,15 @@ import numpy as np from graphqomb.command import TICK, E, M, N -from graphqomb.common import MeasBasis, Plane +from graphqomb.common import Axis, MeasBasis, Plane from graphqomb.gates import MultiGate, SingleGate, TwoQubitGate from graphqomb.pattern import is_runnable from graphqomb.rng import ensure_rng from graphqomb.statevec import StateVector if TYPE_CHECKING: + from numpy.typing import NDArray + from graphqomb.circuit import BaseCircuit from graphqomb.command import Command from graphqomb.gates import Gate @@ -31,6 +33,11 @@ _X_MATRIX = np.asarray([[0, 1], [1, 0]], dtype=np.complex128) _Z_MATRIX = np.asarray([[1, 0], [0, -1]], dtype=np.complex128) +_INPUT_STATE_VECTORS: dict[Axis, NDArray[np.complex128]] = { + Axis.X: np.asarray([1.0, 1.0], dtype=np.complex128) / np.sqrt(2), + Axis.Y: np.asarray([1.0, 1.0j], dtype=np.complex128) / np.sqrt(2), + Axis.Z: np.asarray([1.0, 0.0], dtype=np.complex128), +} class SimulatorBackend(Enum): @@ -146,8 +153,11 @@ def __init__( is_runnable(self.__pattern) if backend == SimulatorBackend.StateVector: - # Note: deterministic check skipped for now - self.state = StateVector.from_num_qubits(len(self.__pattern.input_node_indices)) + input_states = [ + _INPUT_STATE_VECTORS[self.__pattern.input_initialization_axes.get(node, Axis.X)] + for node in self.node_indices + ] + self.state = StateVector.from_product_states(input_states) elif backend == SimulatorBackend.DensityMatrix: raise NotImplementedError else: diff --git a/graphqomb/statevec.py b/graphqomb/statevec.py index 1af2333d..71583cf7 100644 --- a/graphqomb/statevec.py +++ b/graphqomb/statevec.py @@ -121,6 +121,34 @@ def from_num_qubits(num_qubits: int) -> StateVector: raise ValueError(msg) return StateVector(np.full((2,) * num_qubits, 1 / math.sqrt(2**num_qubits), dtype=np.complex128)) + @staticmethod + def from_product_states(states: Sequence[ArrayLike]) -> StateVector: + r"""Create a state vector from one-qubit product states. + + Parameters + ---------- + states : sequence of array-like + One-qubit state vectors in external qubit order. + + Returns + ------- + `StateVector` + The resulting product state vector. + + Raises + ------ + ValueError + If any one-qubit state does not have exactly two amplitudes. + """ + result = np.asarray(1.0, dtype=np.complex128) + for state in states: + one_qubit_state = np.asarray(state, dtype=np.complex128) + if one_qubit_state.size != 2: # noqa: PLR2004 + msg = "Each product-state factor must have exactly two amplitudes." + raise ValueError(msg) + result = np.asarray(np.kron(result, one_qubit_state.reshape(2)), dtype=np.complex128) + return StateVector(result) + @staticmethod def tensor_product(a: StateVector, b: StateVector) -> StateVector: """Tensor product with other state vector, self ⊗ other. diff --git a/graphqomb/stim_compiler.py b/graphqomb/stim_compiler.py index db62cd84..4cc8e16f 100644 --- a/graphqomb/stim_compiler.py +++ b/graphqomb/stim_compiler.py @@ -79,7 +79,8 @@ def _emit_input_nodes(self) -> None: coordinates = self._pattern.input_coordinates if self._emit_qubit_coords else None for node in self._pattern.input_node_indices: coord = coordinates.get(node) if coordinates else None - self._process_prepare(node, coord, is_input=True) + axis = self._pattern.input_initialization_axes.get(node, Axis.X) + self._process_prepare(node, coord, is_input=True, init_axis=axis) def _process_commands(self) -> None: for cmd in self._pattern: @@ -95,7 +96,14 @@ def _process_commands(self) -> None: msg = f"Unsupported command for stim compilation: {type(cmd).__name__}" raise TypeError(msg) - def _process_prepare(self, node: int, coordinate: tuple[float, ...] | None, *, is_input: bool) -> None: + def _process_prepare( + self, + node: int, + coordinate: tuple[float, ...] | None, + *, + is_input: bool, + init_axis: Axis = Axis.X, + ) -> None: event = PrepareEvent(time=self._tick, node=self._node_info(node), is_input=is_input) ops = self._validate_ops_for_event(event, self._collect_noise_ops_from_models(lambda m: m.on_prepare(event))) default_placement = default_noise_placement(event) @@ -104,7 +112,8 @@ def _process_prepare(self, node: int, coordinate: tuple[float, ...] | None, *, i coord = coordinate if self._emit_qubit_coords else None if coord is not None: self._stim_io.write(f"QUBIT_COORDS({', '.join(str(c) for c in coord)}) {node}\n") - self._stim_io.write(f"RX {node}\n") + reset_instr = {Axis.X: "RX", Axis.Y: "RY", Axis.Z: "R"}[init_axis] + self._stim_io.write(f"{reset_instr} {node}\n") self._rec_index += self._emit_noise_ops(ops, NoisePlacement.AFTER, default_placement) self._alive_nodes.add(node) diff --git a/tests/test_graph_compose.py b/tests/test_graph_compose.py index 510a8e00..d0a7e5ac 100644 --- a/tests/test_graph_compose.py +++ b/tests/test_graph_compose.py @@ -4,7 +4,7 @@ import pytest -from graphqomb.common import Plane, PlannerMeasBasis +from graphqomb.common import Axis, Plane, PlannerMeasBasis from graphqomb.graphstate import BaseGraphState, GraphState, compose @@ -170,6 +170,36 @@ def test_compose_preserves_measurement_bases() -> None: assert composed.meas_bases[mapped_node2] == meas_basis +def test_compose_preserves_surviving_input_initialization_axes() -> None: + """Composition preserves initialization axes for inputs that remain inputs.""" + graph1 = GraphState() + g1_in = graph1.add_node() + g1_out = graph1.add_node() + graph1.add_edge(g1_in, g1_out) + graph1.register_input(g1_in, 0, init_axis=Axis.Y) + graph1.register_output(g1_out, 1) + graph1.assign_meas_basis(g1_in, PlannerMeasBasis(Plane.XY, 0.0)) + + graph2 = GraphState() + g2_in_connected = graph2.add_node() + g2_in_survives = graph2.add_node() + g2_out = graph2.add_node() + graph2.add_edge(g2_in_connected, g2_out) + graph2.add_edge(g2_in_survives, g2_out) + graph2.register_input(g2_in_connected, 1, init_axis=Axis.Z) + graph2.register_input(g2_in_survives, 2, init_axis=Axis.Z) + graph2.register_output(g2_out, 3) + graph2.assign_meas_basis(g2_in_connected, PlannerMeasBasis(Plane.XY, 0.0)) + graph2.assign_meas_basis(g2_in_survives, PlannerMeasBasis(Plane.XY, 0.0)) + + composed, node_map1, node_map2 = compose(graph1, graph2) + + assert composed.input_initialization_axes == { + node_map1[g1_in]: Axis.Y, + node_map2[g2_in_survives]: Axis.Z, + } + + def test_compose_full_connection() -> None: """Test compose where all outputs of graph1 connect to inputs of graph2.""" # Create graph1: input [0] -> output [1, 2] diff --git a/tests/test_graphstate.py b/tests/test_graphstate.py index acd80f72..f707e59f 100644 --- a/tests/test_graphstate.py +++ b/tests/test_graphstate.py @@ -7,7 +7,7 @@ import numpy as np import pytest -from graphqomb.common import Plane, PlannerMeasBasis +from graphqomb.common import Axis, Plane, PlannerMeasBasis from graphqomb.euler import LocalClifford from graphqomb.graphstate import GraphState, bipartite_edges, compose, odd_neighbors @@ -102,6 +102,42 @@ def test_add_node_input_output(graph: GraphState) -> None: assert graph.output_node_indices[node_index] == q_index +def test_register_input_defaults_to_x_initialization_axis(graph: GraphState) -> None: + """Input nodes default to positive X-basis initialization.""" + node_index = graph.add_node() + + graph.register_input(node_index, 0) + + assert graph.input_initialization_axes == {node_index: Axis.X} + + +def test_register_input_accepts_pauli_initialization_axis(graph: GraphState) -> None: + """Input nodes can be initialized in a positive Pauli eigenstate.""" + node_index = graph.add_node() + + graph.register_input(node_index, 0, init_axis=Axis.Y) + + assert graph.input_initialization_axes == {node_index: Axis.Y} + + +def test_assign_input_initialization_axis_updates_input(graph: GraphState) -> None: + """Input initialization axes can be assigned after input registration.""" + node_index = graph.add_node() + graph.register_input(node_index, 0) + + graph.assign_input_initialization_axis(node_index, Axis.Z) + + assert graph.input_initialization_axes == {node_index: Axis.Z} + + +def test_assign_input_initialization_axis_rejects_non_input(graph: GraphState) -> None: + """Only registered input nodes can receive an input initialization axis.""" + node_index = graph.add_node() + + with pytest.raises(ValueError, match="Node is not registered as an input node"): + graph.assign_input_initialization_axis(node_index, Axis.Z) + + def test_ensure_node_exists_raises(graph: GraphState) -> None: """Test ensuring a node exists in the graph.""" with pytest.raises(ValueError, match="Node does not exist node=1"): diff --git a/tests/test_graphstate_bulk_init.py b/tests/test_graphstate_bulk_init.py index 17631110..1902dcc7 100644 --- a/tests/test_graphstate_bulk_init.py +++ b/tests/test_graphstate_bulk_init.py @@ -6,7 +6,7 @@ import pytest -from graphqomb.common import Plane, PlannerMeasBasis +from graphqomb.common import Axis, Plane, PlannerMeasBasis from graphqomb.graphstate import GraphState @@ -76,6 +76,22 @@ def test_from_graph_with_multiple_inputs_outputs() -> None: assert gs.output_node_indices == {node_map[outputs[0]]: 0, node_map[outputs[1]]: 1} +def test_from_graph_with_input_initialization_axes() -> None: + """Test from_graph() preserves specified input initialization axes.""" + nodes = ["a", "b", "c"] + edges = [("a", "c"), ("b", "c")] + inputs = ["a", "b"] + + gs, node_map = GraphState.from_graph( + nodes=nodes, + edges=edges, + inputs=inputs, + input_initialization_axes={"a": Axis.Y, "b": Axis.Z}, + ) + + assert gs.input_initialization_axes == {node_map["a"]: Axis.Y, node_map["b"]: Axis.Z} + + def test_from_graph_with_meas_bases() -> None: """Test from_graph() with measurement bases.""" nodes = ["a", "b", "c"] @@ -213,6 +229,19 @@ def test_from_base_graph_state_preserves_indices() -> None: assert dst.output_node_indices[node_map[n1]] == 10 +def test_from_base_graph_state_preserves_input_initialization_axes() -> None: + """Test from_base_graph_state() preserves input initialization axes.""" + src = GraphState() + n0 = src.add_node() + n1 = src.add_node() + src.register_input(n0, 0, init_axis=Axis.Y) + src.register_input(n1, 1, init_axis=Axis.Z) + + dst, node_map = GraphState.from_base_graph_state(src) + + assert dst.input_initialization_axes == {node_map[n0]: Axis.Y, node_map[n1]: Axis.Z} + + def test_from_base_graph_state_independence() -> None: """Test that modifications to copied graph don't affect original.""" # Create source graph diff --git a/tests/test_ptn_format.py b/tests/test_ptn_format.py index ad14ec6d..60521566 100644 --- a/tests/test_ptn_format.py +++ b/tests/test_ptn_format.py @@ -104,6 +104,7 @@ def assert_pattern_equivalent(actual: Pattern, expected: Pattern) -> None: assert actual.input_node_indices == expected.input_node_indices assert actual.output_node_indices == expected.output_node_indices assert actual.input_coordinates == expected.input_coordinates + assert actual.input_initialization_axes == expected.input_initialization_axes assert actual.pauli_frame.xflow == expected.pauli_frame.xflow assert actual.pauli_frame.zflow == expected.pauli_frame.zflow assert actual.pauli_frame.parity_check_group == expected.pauli_frame.parity_check_group @@ -116,13 +117,102 @@ def test_dumps_basic() -> None: pattern = create_simple_pattern() ptn_str = dumps(pattern) - assert ".version 1" in ptn_str + assert ".version 2" in ptn_str assert ".input" in ptn_str assert ".output" in ptn_str assert "#======== QUANTUM ========" in ptn_str assert "#======== CLASSICAL ========" in ptn_str +def test_dumps_omits_default_input_basis() -> None: + """Default X-basis input initialization is implicit in .ptn output.""" + pattern = create_simple_pattern() + + assert ".input_basis" not in dumps(pattern) + + +def test_ptn_roundtrip_preserves_input_initialization_axes() -> None: + """Non-default input initialization axes survive .ptn roundtrip.""" + graph = GraphState() + in0 = graph.add_node() + in1 = graph.add_node() + out0 = graph.add_node() + out1 = graph.add_node() + + graph.register_input(in0, 0, init_axis=Axis.Y) + graph.register_input(in1, 1, init_axis=Axis.Z) + graph.register_output(out0, 0) + graph.register_output(out1, 1) + graph.add_edge(in0, out0) + graph.add_edge(in1, out1) + graph.assign_meas_basis(in0, AxisMeasBasis(Axis.X, Sign.PLUS)) + graph.assign_meas_basis(in1, AxisMeasBasis(Axis.X, Sign.PLUS)) + + pattern = qompile(graph, {in0: {out0}, in1: {out1}}) + ptn_str = dumps(pattern) + loaded = loads(ptn_str) + + assert ".input_basis" in ptn_str + assert loaded.input_initialization_axes == {in0: Axis.Y, in1: Axis.Z} + + +def test_ptn_loads_legacy_input_without_input_basis_as_x() -> None: + """Existing .ptn files without .input_basis default inputs to X initialization.""" + ptn_str = """# GraphQOMB Pattern Format v1 +.version 1 +.input 0:0 +.output 1:0 + +[0] +E 0 1 +M 0 X + + +.xflow 0 -> 1 +""" + + loaded = loads(ptn_str) + + assert loaded.input_initialization_axes == {0: Axis.X} + + +def test_ptn_load_rejects_input_basis_for_non_input_node() -> None: + """Input basis directives can only reference input nodes.""" + ptn_str = """# GraphQOMB Pattern Format v2 +.version 2 +.input 0:0 +.input_basis 1:Z +.output 1:0 + +[0] +E 0 1 +M 0 X + + +.xflow 0 -> 1 +""" + + with pytest.raises(ValueError, match="Input basis specified for non-input node"): + loads(ptn_str) + + +def test_ptn_load_rejects_input_basis_in_legacy_version() -> None: + """Version 1 .ptn files cannot use the version 2 input-basis directive.""" + ptn_str = """# GraphQOMB Pattern Format v1 +.version 1 +.input 0:0 +.input_basis 0:Z +.output 1:0 + +[0] +E 0 1 +M 0 X + + +.xflow 0 -> 1 +""" + + with pytest.raises(ValueError, match=r"\.input_basis requires \.ptn version 2"): + loads(ptn_str) + + def test_dumps_contains_commands() -> None: """Test that dumps includes all command types.""" pattern = create_simple_pattern() diff --git a/tests/test_simulator.py b/tests/test_simulator.py index ec12f2fc..3ae6d7f5 100644 --- a/tests/test_simulator.py +++ b/tests/test_simulator.py @@ -13,10 +13,15 @@ from graphqomb.statevec import StateVector -def _single_output_pattern(*, measured: bool, axis: Axis = Axis.Z) -> tuple[Pattern, int]: +def _single_output_pattern( + *, + measured: bool, + axis: Axis = Axis.Z, + init_axis: Axis = Axis.X, +) -> tuple[Pattern, int]: graph = GraphState() node = graph.add_node() - graph.register_input(node, 0) + graph.register_input(node, 0, init_axis=init_axis) graph.register_output(node, 0) pauli_frame = PauliFrame(graph, xflow={}, zflow={}) @@ -26,10 +31,41 @@ def _single_output_pattern(*, measured: bool, axis: Axis = Axis.Z) -> tuple[Patt output_node_indices=graph.output_node_indices, commands=commands, pauli_frame=pauli_frame, + input_initialization_axes=graph.input_initialization_axes, ) return pattern, node +def test_pattern_simulator_initializes_input_in_x_basis() -> None: + """PatternSimulator initializes X-axis inputs as |+>.""" + pattern, _ = _single_output_pattern(measured=False, init_axis=Axis.X) + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + + simulator.simulate() + + np.testing.assert_allclose(simulator.state.state(), np.asarray([1.0, 1.0]) / np.sqrt(2)) + + +def test_pattern_simulator_initializes_input_in_y_basis() -> None: + """PatternSimulator initializes Y-axis inputs as |Y+>.""" + pattern, _ = _single_output_pattern(measured=False, init_axis=Axis.Y) + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + + simulator.simulate() + + np.testing.assert_allclose(simulator.state.state(), np.asarray([1.0, 1.0j]) / np.sqrt(2)) + + +def test_pattern_simulator_initializes_input_in_z_basis() -> None: + """PatternSimulator initializes Z-axis inputs as |0>.""" + pattern, _ = _single_output_pattern(measured=False, init_axis=Axis.Z) + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + + simulator.simulate() + + np.testing.assert_allclose(simulator.state.state(), np.asarray([1.0, 0.0])) + + def test_pattern_simulator_applies_output_x_frame_to_statevector() -> None: """An unmeasured output statevector should include pending X frame corrections.""" pattern, node = _single_output_pattern(measured=False) diff --git a/tests/test_stim_compiler.py b/tests/test_stim_compiler.py index 116c76e3..40ddfc33 100644 --- a/tests/test_stim_compiler.py +++ b/tests/test_stim_compiler.py @@ -136,6 +136,52 @@ def test_stim_compile_basic_pattern() -> None: assert stim_str.count("\n") > 0 +@pytest.mark.parametrize( + ("init_axis", "expected_reset"), + [ + (Axis.X, "RX"), + (Axis.Y, "RY"), + (Axis.Z, "R"), + ], +) +def test_stim_compile_uses_input_initialization_axis(init_axis: Axis, expected_reset: str) -> None: + """Input initialization axes choose the corresponding Stim reset instruction.""" + graph = GraphState() + in_node = graph.add_node() + out_node = graph.add_node() + + graph.register_input(in_node, 0, init_axis=init_axis) + graph.register_output(out_node, 0) + graph.add_edge(in_node, out_node) + graph.assign_meas_basis(in_node, AxisMeasBasis(Axis.X, Sign.PLUS)) + + pattern = qompile(graph, {in_node: {out_node}}) + stim_lines = stim_compile(pattern).splitlines() + + assert f"{expected_reset} {in_node}" in stim_lines + + +def test_stim_compile_keeps_non_input_preparations_in_x_basis() -> None: + """Only input reset instructions use the input initialization axis.""" + graph = GraphState() + in_node = graph.add_node() + mid_node = graph.add_node() + out_node = graph.add_node() + + graph.register_input(in_node, 0, init_axis=Axis.Z) + graph.register_output(out_node, 0) + graph.add_edge(in_node, mid_node) + graph.add_edge(mid_node, out_node) + graph.assign_meas_basis(in_node, AxisMeasBasis(Axis.X, Sign.PLUS)) + graph.assign_meas_basis(mid_node, AxisMeasBasis(Axis.X, Sign.PLUS)) + + pattern = qompile(graph, {in_node: {mid_node}, mid_node: {out_node}}) + stim_lines = stim_compile(pattern).splitlines() + + assert f"R {in_node}" in stim_lines + assert f"RX {mid_node}" in stim_lines + + def test_stim_compile_x_measurement() -> None: """Test X measurement compilation.""" pattern, meas_node, in_node = create_simple_pattern_x_measurement()