diff --git a/CHANGELOG.md b/CHANGELOG.md index 44cabffa..5470441c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ 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. +- **Stim Circuit Import**: Added `stim_file_to_pattern()`, `stim_text_to_pattern()`, and `stim_circuit_to_pattern()` for converting supported Stim circuits into GraphQOMB patterns. The importer handles Clifford unitary blocks and `TICK`-separated Pauli measurements (`M`/`MZ`, `MX`, `MY`, `MXX`, `MYY`, `MZZ`, and `MPP`), assigns single-qubit measurement bases directly to their graph nodes, terminates each qubit lifetime at its direct measurement while allowing disjoint qubits to continue, validates that same-block MPP products commute, supports Type I and Type II Y foliation, causally lowers commuting MPP groups, composes each MPP unit through a separate unmeasured output layer, and resolves detector/logical-observable records once across the full flattened circuit. Circuit-level noise and measurement-error probabilities are intentionally omitted because GraphQOMB uses an MBQC-specific noise model; reset, measurement-reset, and feedback instructions remain unsupported. ### Changed @@ -17,6 +18,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - **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. +### Fixed + +- **Stim Compiler**: Preserve minus-signed axis measurement bases using inverted Stim measurement targets. +- **State Vector Array Conversion**: Convert to a real NumPy dtype without warnings when every amplitude is real-valued, and reject the conversion when it would discard nonzero imaginary amplitudes. + ### Removed - **Pattern Commands**: Removed explicit `graphqomb.command.X` and `graphqomb.command.Z` correction commands. Output corrections are now represented by `PauliFrame` only, and `.ptn` files containing `X`/`Z` command lines are rejected. diff --git a/README.md b/README.md index dce030ce..07824009 100644 --- a/README.md +++ b/README.md @@ -115,6 +115,7 @@ If you already have a graph-state design and explicit feedforward maps, you can - **API reference**: https://graphqomb.readthedocs.io/en/latest/references.html - **QEC graph-state builder reference**: https://graphqomb.readthedocs.io/en/latest/qec.html - **Stim MPP import reference**: https://graphqomb.readthedocs.io/en/latest/stim_mpp.html +- **Stim circuit import reference**: https://graphqomb.readthedocs.io/en/latest/stim_importer.html - **Stim compiler reference**: https://graphqomb.readthedocs.io/en/latest/stim_compiler.html ## Current Scope diff --git a/docs/source/qec.rst b/docs/source/qec.rst index e0c33f3d..45052d01 100644 --- a/docs/source/qec.rst +++ b/docs/source/qec.rst @@ -9,6 +9,13 @@ graph-state builder uses the first two qubit-coordinate components for the data plane and the first three explicitly supplied stabilizer-coordinate components for ancilla placement. +With ``data_as_io=True``, the stabilizer-measurement unit has a separate, +unmeasured output layer. Type I therefore has two measured data layers followed +by an output layer. Type II uses three Y-measured layers for qubits with Y +support, followed by an output layer; qubits without Y support retain the Type I +two-measurement-layer layout. Ancilla support edges only touch measurement +layers, and the output of one composed unit becomes the input of the next. + .. automodule:: graphqomb.qec.qeccode :members: :show-inheritance: diff --git a/docs/source/references.rst b/docs/source/references.rst index 1d252c18..05fc9ae0 100644 --- a/docs/source/references.rst +++ b/docs/source/references.rst @@ -23,6 +23,7 @@ Module reference scheduler qec stim_mpp + stim_importer stim_compiler noise_model visualizer diff --git a/docs/source/stim_importer.rst b/docs/source/stim_importer.rst new file mode 100644 index 00000000..479ad4fa --- /dev/null +++ b/docs/source/stim_importer.rst @@ -0,0 +1,84 @@ +Stim circuit import +=================== + +Install the optional Stim integration before importing this module: + +.. code-block:: console + + uv add "graphqomb[stim]" + +The circuit importer converts supported Stim circuits into GraphQOMB +measurement patterns. It accepts Clifford unitary blocks and Pauli measurement +blocks separated by ``TICK``. The supported Pauli measurement instructions are +``M``/``MZ``, ``MX``, ``MY``, ``MXX``, ``MYY``, ``MZZ``, and ``MPP``. + +Single-qubit measurements assign an ``AxisMeasBasis`` directly to the measured +graph node. They do not create an ``MPP`` extraction or an ancillary parity +measurement node. Inverted single-qubit measurement targets select the minus +sign of that node's basis. A direct single-qubit measurement terminates that +qubit's lifetime: a later quantum operation on the same qubit is rejected, +while operations on other qubits may continue. Reset and qubit reuse are not +supported. + +Two-qubit measurements are parity measurements and are lowered to equivalent +unsigned ``MPP`` products. Inverted targets in ``MXX``, ``MYY``, ``MZZ``, and +``MPP`` are rejected because GraphQOMB does not currently retain the +corresponding parity offset. + +All ``MPP`` instructions within one ``TICK`` block are represented by one +combined extraction and are validated to commute. Anticommuting products in +the same block are rejected. The importer uses one compact +stabilizer-measurement unit when its correction flow is causal. If the compact +unit has a cyclic Pauli flow, the commuting products are lowered to equivalent +sequential units instead. Non-commuting measurements must be separated by +``TICK`` in the source circuit. Each unit has a distinct unmeasured output +layer, which is composed with the next unitary or measurement fragment by qubit +index. Pass ``y_foliation=YFoliation.TYPE_II`` to any of the three import entry +points to use the three-layer Y-measurement construction; Type I is the default. + +The flattened ideal circuit is analyzed once for measurement records. Records +from single-qubit measurements, pair measurements, ``MPP``, and ideal-zero +``MPAD`` replacements share one absolute index space. ``DETECTOR`` and +``OBSERVABLE_INCLUDE`` targets are resolved against that whole-circuit index +space, and each ``StimMppExtraction`` retains the corresponding absolute record +sets even when an annotation also references records outside that MPP unit. + +Noise policy +------------ + +Circuit-level Stim noise is intentionally not imported. GraphQOMB applies +noise to the compiled MBQC pattern through its own noise-model API, where +preparation, entanglement, measurement, and idle events differ from the source +circuit operations. + +Pure noise instructions are omitted during import. Error probabilities attached +to Pauli measurements are also omitted while retaining the ideal measurement. +Heralded noise records are retained as ideal zero-valued record positions so +that later ``DETECTOR`` and ``OBSERVABLE_INCLUDE`` references remain aligned. + +Reset instructions (``R``/``RZ``, ``RX``, and ``RY``) and combined +measurement-reset instructions (``MR``/``MRZ``, ``MRX``, and ``MRY``) are not +handled by this importer. Consequently, a directly measured qubit cannot begin +a new lifetime later in the circuit. + +.. code-block:: python + + from graphqomb.qec.qeccode import YFoliation + from graphqomb.stim_importer import stim_text_to_pattern + + result = stim_text_to_pattern( + """ + MX 0 + MYY 1 2 + DETECTOR rec[-2] rec[-1] + """, + y_foliation=YFoliation.TYPE_II, + ) + pattern = result.pattern + +API reference +------------- + +.. automodule:: graphqomb.stim_importer + :members: + :show-inheritance: diff --git a/docs/source/stim_mpp.rst b/docs/source/stim_mpp.rst index e10c538e..4e50c5dc 100644 --- a/docs/source/stim_mpp.rst +++ b/docs/source/stim_mpp.rst @@ -19,8 +19,9 @@ ancilla nodes. Set ``coord_dims`` to the number of leading ``QUBIT_COORDS`` components to retain. The importer does not restrict this value to two or three dimensions; -each referenced coordinate declaration must contain at least the requested -number of components. +each referenced qubit's final accumulated coordinate must contain at least the +requested number of components. Repeated declarations are combined according +to Stim's coordinate semantics. .. code-block:: python @@ -34,6 +35,10 @@ number of components. API reference ------------- +.. autoclass:: graphqomb.qec.stim_mpp.StimMppExtraction + :members: + :show-inheritance: + .. automodule:: graphqomb.qec.stim_mpp :members: :show-inheritance: diff --git a/graphqomb/qec/_stim.py b/graphqomb/qec/_stim.py new file mode 100644 index 00000000..1259368f --- /dev/null +++ b/graphqomb/qec/_stim.py @@ -0,0 +1,342 @@ +"""Shared internal helpers for parsing Stim QEC data.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +import numpy as np +from scipy.sparse import csr_array, lil_array + +from graphqomb.qec.qeccode import Coordinate, StabilizerCode + +if TYPE_CHECKING: + from collections.abc import Mapping, Sequence + + import stim + + +PauliSupport = tuple[tuple[int, str], ...] + + +@dataclass(frozen=True) +class StimMppExtraction: + """Stabilizer-code data extracted from Stim MPP products. + + Attributes + ---------- + code : StabilizerCode + Dense-column stabilizer code using the ``[Hx | Hz]`` convention. + stim_to_column : dict[int, int] + Mapping from original Stim qubit ids to dense matrix columns. + column_to_stim : dict[int, int] + Inverse dense-column mapping. + supports : tuple[PauliSupport, ...] + Original Stim Pauli supports, one support per stabilizer row. + detector_rows : tuple[frozenset[int], ...] + Detector groups as selected-MPP stabilizer row indices. + logical_observable_rows : dict[int, frozenset[int]] + Logical observables as selected-MPP stabilizer row indices. + detector_record_indices : tuple[frozenset[int], ...] + Absolute Stim measurement-record indices for selected detectors. + logical_observable_record_indices : dict[int, frozenset[int]] + Absolute Stim record indices for selected logical observables. + """ + + code: StabilizerCode + stim_to_column: dict[int, int] + column_to_stim: dict[int, int] + supports: tuple[PauliSupport, ...] + detector_rows: tuple[frozenset[int], ...] + logical_observable_rows: dict[int, frozenset[int]] + detector_record_indices: tuple[frozenset[int], ...] = () + logical_observable_record_indices: dict[int, frozenset[int]] = field(default_factory=dict) + + def detector_groups(self, ancilla_nodes: Mapping[int, int]) -> list[set[int]]: + """Return detector groups mapped to graph node ids for ``qompile``. + + Returns + ------- + list[set[int]] + Detector groups suitable for ``qompile``. + """ + return [_map_rows_to_nodes(rows, ancilla_nodes, "detector") for rows in self.detector_rows] + + def logical_observables(self, ancilla_nodes: Mapping[int, int]) -> dict[int, set[int]]: + """Return logical observables mapped to graph node ids for ``qompile``. + + Returns + ------- + dict[int, set[int]] + Logical-observable node groups keyed by Stim observable index. + """ + return { + logical_idx: _map_rows_to_nodes(rows, ancilla_nodes, f"logical observable {logical_idx}") + for logical_idx, rows in self.logical_observable_rows.items() + } + + +def stim_mpp_extraction_from_records( + supports: Sequence[PauliSupport], + record_indices: Sequence[int], + *, + coordinate_by_stim_id: Mapping[int, Coordinate], + detector_record_indices: Sequence[frozenset[int]], + logical_observable_record_indices: Mapping[int, frozenset[int]], +) -> StimMppExtraction: + """Build an MPP extraction from globally indexed measurement records. + + Returns + ------- + StimMppExtraction + Extracted stabilizer data and record metadata. + + Raises + ------ + ValueError + If the support and record counts differ. + """ + if len(supports) != len(record_indices): + msg = "MPP support count does not match its measurement-record count." + raise ValueError(msg) + + record_to_row = {record_index: row for row, record_index in enumerate(record_indices)} + selected_detector_rows: list[frozenset[int]] = [] + selected_detector_records: list[frozenset[int]] = [] + for records in detector_record_indices: + rows = frozenset(record_to_row[record] for record in records if record in record_to_row) + if rows: + selected_detector_rows.append(rows) + selected_detector_records.append(records) + + selected_logical_rows: dict[int, frozenset[int]] = {} + selected_logical_records: dict[int, frozenset[int]] = {} + for logical_idx, records in logical_observable_record_indices.items(): + rows = frozenset(record_to_row[record] for record in records if record in record_to_row) + if rows: + selected_logical_rows[logical_idx] = rows + selected_logical_records[logical_idx] = records + + matrix, stim_to_column, column_to_stim, qubit_coords = _build_stabilizer_data( + supports, + coordinate_by_stim_id, + ) + + return StimMppExtraction( + code=StabilizerCode(matrix, qubit_coords=qubit_coords), + stim_to_column=stim_to_column, + column_to_stim=column_to_stim, + supports=tuple(supports), + detector_rows=tuple(selected_detector_rows), + logical_observable_rows=selected_logical_rows, + detector_record_indices=tuple(selected_detector_records), + logical_observable_record_indices=selected_logical_records, + ) + + +def extract_qubit_coordinates( + circuit: stim.Circuit, + *, + coord_dims: int, +) -> dict[int, Coordinate]: + """Return final Stim qubit coordinates projected to ``coord_dims``. + + Returns + ------- + dict[int, Coordinate] + Final coordinates keyed by Stim qubit id. + + Raises + ------ + ValueError + If a coordinate has fewer dimensions than requested. + """ + coordinates: dict[int, Coordinate] = {} + for stim_id, values in circuit.get_final_qubit_coordinates().items(): + if len(values) < coord_dims: + msg = ( + f"QUBIT_COORDS for qubit {stim_id} has {len(values)} coordinate(s), " + f"fewer than requested coord_dims={coord_dims}." + ) + raise ValueError(msg) + coordinates[int(stim_id)] = tuple(float(value) for value in values[:coord_dims]) + return coordinates + + +def record_targets_to_absolute_indices( + targets: Sequence[stim.GateTarget], + *, + measurement_count: int, + instruction_name: str, +) -> frozenset[int]: + """Resolve relative Stim record targets to absolute parity indices. + + Returns + ------- + frozenset[int] + Absolute record indices after parity cancellation. + + Raises + ------ + ValueError + If a target is invalid or references a record before time began. + """ + record_indices: set[int] = set() + for target in targets: + if not target.is_measurement_record_target: + msg = f"{instruction_name} contains unsupported target {target!r}; only rec targets are supported." + raise ValueError(msg) + record_index = measurement_count + int(target.value) + if not 0 <= record_index < measurement_count: + msg = f"{instruction_name} refers to measurement record {record_index} before the beginning of time." + raise ValueError(msg) + if record_index in record_indices: + record_indices.remove(record_index) + else: + record_indices.add(record_index) + return frozenset(record_indices) + + +def observable_index(instruction: stim.CircuitInstruction) -> int: + """Return the logical-observable index from a Stim annotation. + + Returns + ------- + int + Logical-observable index. + + Raises + ------ + ValueError + If the annotation does not have one integer argument. + """ + args = instruction.gate_args_copy() + if len(args) != 1 or not args[0].is_integer(): + msg = "OBSERVABLE_INCLUDE must have one integer observable index." + raise ValueError(msg) + return int(args[0]) + + +def mpp_targets_to_products(targets: Sequence[stim.GateTarget]) -> list[PauliSupport]: + """Parse Stim MPP targets into unsigned Pauli products. + + Returns + ------- + list[PauliSupport] + Parsed Pauli products in target order. + + Raises + ------ + ValueError + If the target sequence is signed or malformed. + """ + products: list[PauliSupport] = [] + current: list[tuple[int, str]] = [] + seen_in_current: set[int] = set() + expect_pauli = True + + for target in targets: + if target.is_combiner: + if expect_pauli: + msg = "Invalid MPP target list: unexpected combiner." + raise ValueError(msg) + expect_pauli = True + continue + + if target.is_inverted_result_target: + msg = "Signed MPP products are not supported; inverted Pauli targets cannot be imported." + raise ValueError(msg) + pauli = _target_pauli(target) + if current and not expect_pauli: + products.append(tuple(current)) + current = [] + seen_in_current = set() + + qid = int(target.value) + if qid in seen_in_current: + msg = f"Invalid MPP product: qubit {qid} appears more than once." + raise ValueError(msg) + current.append((qid, pauli)) + seen_in_current.add(qid) + expect_pauli = False + + if expect_pauli: + msg = "Invalid MPP target list: trailing combiner or empty product." + raise ValueError(msg) + products.append(tuple(current)) + return products + + +def pauli_products_commute(left: PauliSupport, right: PauliSupport) -> bool: + """Return whether two unsigned Pauli products commute. + + Returns + ------- + bool + Whether the two products commute. + """ + right_by_qubit = dict(right) + anticommuting_overlaps = sum(qubit in right_by_qubit and pauli != right_by_qubit[qubit] for qubit, pauli in left) + return anticommuting_overlaps % 2 == 0 + + +def plain_qubit_target(target: stim.GateTarget, instruction_name: str) -> int: + """Return a plain Stim qubit target. + + Returns + ------- + int + Stim qubit id. + + Raises + ------ + ValueError + If the target is not a plain qubit target. + """ + qubit_value = target.qubit_value + if qubit_value is None or not target.is_qubit_target: + msg = f"{instruction_name} contains unsupported target {target!r}; only plain qubit targets are supported." + raise ValueError(msg) + return int(qubit_value) + + +def _build_stabilizer_data( + supports: Sequence[PauliSupport], + coordinate_by_stim_id: Mapping[int, Coordinate], +) -> tuple[csr_array[Any, tuple[int, int]], dict[int, int], dict[int, int], dict[int, Coordinate]]: + stim_ids = sorted({qid for support in supports for qid, _pauli in support}) + stim_to_column = {qid: column for column, qid in enumerate(stim_ids)} + column_to_stim = {column: qid for qid, column in stim_to_column.items()} + + num_qubits = len(stim_ids) + matrix = lil_array((len(supports), 2 * num_qubits), dtype=np.bool_) + for row, support in enumerate(supports): + for stim_id, pauli in support: + column = stim_to_column[stim_id] + if pauli in {"X", "Y"}: + matrix[row, column] = True + if pauli in {"Z", "Y"}: + matrix[row, num_qubits + column] = True + + qubit_coords = {stim_to_column[qid]: coord for qid, coord in coordinate_by_stim_id.items() if qid in stim_to_column} + stabilizer_matrix = csr_array(matrix, shape=(len(supports), 2 * num_qubits)) + return stabilizer_matrix, stim_to_column, column_to_stim, qubit_coords + + +def _target_pauli(target: stim.GateTarget) -> str: + if target.is_x_target: + return "X" + if target.is_y_target: + return "Y" + if target.is_z_target: + return "Z" + msg = f"Unsupported MPP target: {target!r}." + raise ValueError(msg) + + +def _map_rows_to_nodes(rows: frozenset[int], ancilla_nodes: Mapping[int, int], label: str) -> set[int]: + missing_rows = sorted(row for row in rows if row not in ancilla_nodes) + if missing_rows: + msg = f"Cannot map {label}; ancilla node map is missing stabilizer row(s): {missing_rows}." + raise ValueError(msg) + return {ancilla_nodes[row] for row in rows} diff --git a/graphqomb/qec/qeccode.py b/graphqomb/qec/qeccode.py index c8c5be86..22c25b73 100644 --- a/graphqomb/qec/qeccode.py +++ b/graphqomb/qec/qeccode.py @@ -71,9 +71,12 @@ class _DataLayerPlan(NamedTuple): """Data-node layer layout for graph-state construction.""" data_layers: dict[int, tuple[int, ...]] + measurement_layers: dict[int, tuple[int, ...]] coordinate_z_by_layer: dict[tuple[int, int], float] meas_basis_by_qubit: dict[int, AxisMeasBasis] y_foliation: YFoliation + data_as_io: bool + qubit_indices: Mapping[int, int] | None class _StabilizerSupport(NamedTuple): @@ -88,6 +91,8 @@ def build_graph_state( z_base: int = 0, *, y_foliation: YFoliation = YFoliation.TYPE_I, + data_as_io: bool = False, + qubit_indices: Mapping[int, int] | None = None, ) -> StabilizerGraphStateBuildResult: """Build a graph-state unit from a stabilizer code. @@ -97,11 +102,18 @@ def build_graph_state( Stabilizer code to convert. The X support is connected to the upper data layer and the Z support is connected to the lower data layer. z_base : `int`, optional - Lower data-layer index. The builder creates layers ``z_base`` and - ``z_base + 1``, by default 0. + Lower stabilizer-measurement data-layer index, by default 0. Type I + uses two measurement layers; Type II uses three for Y support. When + ``data_as_io`` is enabled, a separate output layer is appended. y_foliation : `YFoliation`, optional Foliation variant. Type II uses a three-node Y-measured data chain only for qubits that have an Hx=Hz=1 support in at least one stabilizer row. + data_as_io : `bool`, optional + Whether to register the first stabilizer-measurement data nodes as + inputs and append separate unmeasured output nodes, by default False. + qubit_indices : collections.abc.Mapping[int, int] | None, optional + Mapping from stabilizer-code qubit columns to graph qindices when + ``data_as_io`` is enabled. If omitted, code qubit columns are used. Returns ------- @@ -112,10 +124,30 @@ def build_graph_state( ------ TypeError If z_base is not an integer. + ValueError + If ``qubit_indices`` is invalid for the requested data I/O layout. """ if not isinstance(z_base, int): msg = "z_base must be an integer." raise TypeError(msg) + if qubit_indices is not None: + if not data_as_io: + msg = "qubit_indices can only be used when data_as_io=True." + raise ValueError(msg) + expected_qubits = set(range(code.num_qubits)) + provided_qubits = set(qubit_indices) + if provided_qubits != expected_qubits: + missing = sorted(expected_qubits - provided_qubits) + unexpected = sorted(provided_qubits - expected_qubits) + msg = ( + "qubit_indices must map every stabilizer-code qubit exactly once; " + f"missing={missing}, unexpected={unexpected}." + ) + raise ValueError(msg) + qindices = list(qubit_indices.values()) + if len(qindices) != len(set(qindices)): + msg = "qubit_indices values must be unique." + raise ValueError(msg) graph = GraphState() x_meas_basis = AxisMeasBasis(Axis.X, Sign.PLUS) @@ -124,6 +156,8 @@ def build_graph_state( code, z_base=z_base, y_foliation=y_foliation, + data_as_io=data_as_io, + qubit_indices=qubit_indices, ) data_nodes = _add_layered_data_nodes(graph, code, data_layer_plan) ancilla_nodes = _add_ancilla_nodes(graph, code, data_nodes, data_layer_plan, x_meas_basis) @@ -136,8 +170,11 @@ def _data_layer_plan( *, z_base: int, y_foliation: YFoliation, + data_as_io: bool, + qubit_indices: Mapping[int, int] | None, ) -> _DataLayerPlan: data_layers: dict[int, tuple[int, ...]] = {} + measurement_layers: dict[int, tuple[int, ...]] = {} coordinate_z_by_layer: dict[tuple[int, int], float] = {} meas_basis_by_qubit: dict[int, AxisMeasBasis] = {} x_meas_basis = AxisMeasBasis(Axis.X, Sign.PLUS) @@ -145,31 +182,42 @@ def _data_layer_plan( y_chain_qubits: set[int] = _qubits_with_y_support(code) if y_foliation is YFoliation.TYPE_II else set() for qubit in range(code.num_qubits): - layers: tuple[int, ...] - coordinate_zs: tuple[float, ...] + measured_layers: tuple[int, ...] + measured_coordinate_zs: tuple[float, ...] if qubit in y_chain_qubits: - layers = (z_base, z_base + 1, z_base + 2) - coordinate_zs = ( + measured_layers = (z_base, z_base + 1, z_base + 2) + measured_coordinate_zs = ( float(z_base), float(z_base) + 0.5, float(z_base + 1), ) meas_basis = y_meas_basis else: - layers = (z_base, z_base + 1) - coordinate_zs = (float(z_base), float(z_base + 1)) + measured_layers = (z_base, z_base + 1) + measured_coordinate_zs = (float(z_base), float(z_base + 1)) meas_basis = x_meas_basis + if data_as_io: + layers = (*measured_layers, measured_layers[-1] + 1) + coordinate_zs = (*measured_coordinate_zs, float(z_base + 2)) + else: + layers = measured_layers + coordinate_zs = measured_coordinate_zs + data_layers[qubit] = layers + measurement_layers[qubit] = measured_layers meas_basis_by_qubit[qubit] = meas_basis for layer, coordinate_z in zip(layers, coordinate_zs, strict=True): coordinate_z_by_layer[qubit, layer] = coordinate_z return _DataLayerPlan( data_layers=data_layers, + measurement_layers=measurement_layers, coordinate_z_by_layer=coordinate_z_by_layer, meas_basis_by_qubit=meas_basis_by_qubit, y_foliation=y_foliation, + data_as_io=data_as_io, + qubit_indices=qubit_indices, ) @@ -188,16 +236,24 @@ def _add_layered_data_nodes( data_nodes: dict[tuple[int, int], int] = {} for qubit in range(code.num_qubits): previous_node: int | None = None - for layer in data_layer_plan.data_layers[qubit]: + layers = data_layer_plan.data_layers[qubit] + measured_layers = set(data_layer_plan.measurement_layers[qubit]) + for layer in layers: node = graph.add_node( coordinate=_data_coordinate(code, qubit, data_layer_plan.coordinate_z_by_layer[qubit, layer]) ) if previous_node is not None: graph.add_edge(previous_node, node) - graph.assign_meas_basis(node, data_layer_plan.meas_basis_by_qubit[qubit]) + if layer in measured_layers: + graph.assign_meas_basis(node, data_layer_plan.meas_basis_by_qubit[qubit]) data_nodes[qubit, layer] = node previous_node = node + if data_layer_plan.data_as_io: + q_index = data_layer_plan.qubit_indices[qubit] if data_layer_plan.qubit_indices is not None else qubit + graph.register_input(data_nodes[qubit, data_layer_plan.measurement_layers[qubit][0]], q_index) + graph.register_output(data_nodes[qubit, layers[-1]], q_index) + return data_nodes @@ -255,7 +311,7 @@ def _connect_stabilizer_support( ) -> list[int]: connected_data_nodes: list[int] = [] for qubit in sorted(support.hx | support.hz): - layers = data_layer_plan.data_layers[qubit] + layers = data_layer_plan.measurement_layers[qubit] if data_layer_plan.y_foliation is YFoliation.TYPE_II and len(layers) == _TYPE_II_CHAIN_LENGTH: layer = _type_ii_support_layer(layers, has_x=qubit in support.hx, has_z=qubit in support.hz) data_node = data_nodes[qubit, layer] diff --git a/graphqomb/qec/stim_mpp.py b/graphqomb/qec/stim_mpp.py index 153e6b7d..de000b73 100644 --- a/graphqomb/qec/stim_mpp.py +++ b/graphqomb/qec/stim_mpp.py @@ -2,21 +2,24 @@ from __future__ import annotations -from dataclasses import dataclass, field +from dataclasses import dataclass from pathlib import Path -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING -import numpy as np import stim -from scipy.sparse import csr_array, lil_array -from graphqomb.qec.qeccode import Coordinate, StabilizerCode +from graphqomb.qec._stim import ( + PauliSupport, + StimMppExtraction, + extract_qubit_coordinates, + mpp_targets_to_products, + observable_index, + record_targets_to_absolute_indices, + stim_mpp_extraction_from_records, +) if TYPE_CHECKING: - from collections.abc import Mapping, Sequence - - -PauliSupport = tuple[tuple[int, str], ...] + from collections.abc import Sequence @dataclass(frozen=True) @@ -26,82 +29,7 @@ class _MppProductRecord: @dataclass(frozen=True) -class StimMppExtraction: - """Stabilizer-code data extracted from Stim MPP products. - - Attributes - ---------- - code : StabilizerCode - Dense-column stabilizer code using the ``[Hx | Hz]`` convention. - stim_to_column : dict[int, int] - Mapping from original Stim qubit ids to dense matrix columns. - column_to_stim : dict[int, int] - Inverse dense-column mapping. - supports : tuple[PauliSupport, ...] - Original Stim Pauli supports, one support per stabilizer row. - detector_rows : tuple[frozenset[int], ...] - Detector groups as selected-MPP stabilizer row indices. If a Stim - detector also references measurements outside the selected MPP products, - only rows represented in this extraction are included here. - logical_observable_rows : dict[int, frozenset[int]] - Logical observables as selected-MPP stabilizer row indices, keyed by - Stim logical observable index. External measurement records are ignored - in this row view. - detector_record_indices : tuple[frozenset[int], ...] - Absolute Stim measurement-record indices for detectors that touch at - least one MPP product represented in this extraction. - logical_observable_record_indices : dict[int, frozenset[int]] - Absolute Stim measurement-record indices for logical observables that - touch at least one MPP product represented in this extraction. - """ - - code: StabilizerCode - stim_to_column: dict[int, int] - column_to_stim: dict[int, int] - supports: tuple[PauliSupport, ...] - detector_rows: tuple[frozenset[int], ...] - logical_observable_rows: dict[int, frozenset[int]] - detector_record_indices: tuple[frozenset[int], ...] = () - logical_observable_record_indices: dict[int, frozenset[int]] = field(default_factory=dict) - - def detector_groups(self, ancilla_nodes: Mapping[int, int]) -> list[set[int]]: - """Return detector groups mapped to graph node ids for ``qompile``. - - Parameters - ---------- - ancilla_nodes : collections.abc.Mapping[int, int] - Mapping from selected-MPP stabilizer rows to graph node ids. - - Returns - ------- - list[set[int]] - Detector groups suitable for ``qompile(..., parity_check_group=...)``. - """ - return [_map_rows_to_nodes(rows, ancilla_nodes, "detector") for rows in self.detector_rows] - - def logical_observables(self, ancilla_nodes: Mapping[int, int]) -> dict[int, set[int]]: - """Return logical observables mapped to graph node ids for ``qompile``. - - Parameters - ---------- - ancilla_nodes : collections.abc.Mapping[int, int] - Mapping from selected-MPP stabilizer rows to graph node ids. - - Returns - ------- - dict[int, set[int]] - Logical observables suitable for ``qompile(..., logical_observables=...)``. - """ - return { - logical_idx: _map_rows_to_nodes(rows, ancilla_nodes, f"logical observable {logical_idx}") - for logical_idx, rows in self.logical_observable_rows.items() - } - - -@dataclass(frozen=True) -class _StimMppAnnotations: - detector_rows: tuple[frozenset[int], ...] - logical_observable_rows: dict[int, frozenset[int]] +class _StimRecordAnnotations: detector_record_indices: tuple[frozenset[int], ...] logical_observable_record_indices: dict[int, frozenset[int]] @@ -161,7 +89,7 @@ def stabilizer_code_from_stim_text( msg = "mpp_layer must be non-negative." raise ValueError(msg) circuit = stim.Circuit(text).flattened() - coordinate_by_stim_id = _extract_qubit_coordinates(circuit, coord_dims=coord_dims) + coordinate_by_stim_id = extract_qubit_coordinates(circuit, coord_dims=coord_dims) layers = _extract_mpp_layers(circuit) selected_layer = _select_mpp_products(layers, mpp_layer=mpp_layer) @@ -170,74 +98,16 @@ def stabilizer_code_from_stim_text( layer_label = "file" if mpp_layer is None else f"layer {mpp_layer}" msg = f"MPP {layer_label} is empty." raise ValueError(msg) - record_to_row = {product.record_index: row for row, product in enumerate(selected_layer)} - annotations = _extract_selected_mpp_annotations( - circuit, - record_to_row=record_to_row, - ) - - matrix, stim_to_column, column_to_stim, qubit_coords = _build_stabilizer_data( + annotations = _extract_stim_record_annotations(circuit) + return stim_mpp_extraction_from_records( supports, - coordinate_by_stim_id, - ) - - return StimMppExtraction( - code=StabilizerCode(matrix, qubit_coords=qubit_coords), - stim_to_column=stim_to_column, - column_to_stim=column_to_stim, - supports=supports, - detector_rows=annotations.detector_rows, - logical_observable_rows=annotations.logical_observable_rows, + tuple(product.record_index for product in selected_layer), + coordinate_by_stim_id=coordinate_by_stim_id, detector_record_indices=annotations.detector_record_indices, logical_observable_record_indices=annotations.logical_observable_record_indices, ) -def _build_stabilizer_data( - supports: Sequence[PauliSupport], - coordinate_by_stim_id: Mapping[int, Coordinate], -) -> tuple[csr_array[Any, tuple[int, int]], dict[int, int], dict[int, int], dict[int, Coordinate]]: - stim_ids = sorted({qid for support in supports for qid, _pauli in support}) - stim_to_column = {qid: column for column, qid in enumerate(stim_ids)} - column_to_stim = {column: qid for qid, column in stim_to_column.items()} - - num_qubits = len(stim_ids) - matrix = lil_array((len(supports), 2 * num_qubits), dtype=np.bool_) - for row, support in enumerate(supports): - for stim_id, pauli in support: - column = stim_to_column[stim_id] - if pauli in {"X", "Y"}: - matrix[row, column] = True - if pauli in {"Z", "Y"}: - matrix[row, num_qubits + column] = True - - qubit_coords = {stim_to_column[qid]: coord for qid, coord in coordinate_by_stim_id.items() if qid in stim_to_column} - stabilizer_matrix = csr_array(matrix, shape=(len(supports), 2 * num_qubits)) - return stabilizer_matrix, stim_to_column, column_to_stim, qubit_coords - - -def _extract_qubit_coordinates( - circuit: stim.Circuit, - *, - coord_dims: int, -) -> dict[int, Coordinate]: - coordinates: dict[int, Coordinate] = {} - for instruction in circuit: - if not isinstance(instruction, stim.CircuitInstruction): - msg = "Flattened Stim circuit unexpectedly contains a repeat block." - raise TypeError(msg) - if instruction.name != "QUBIT_COORDS": - continue - args = instruction.gate_args_copy() - if len(args) < coord_dims: - msg = f"QUBIT_COORDS has {len(args)} coordinate(s), fewer than requested coord_dims={coord_dims}." - raise ValueError(msg) - coord = tuple(float(value) for value in args[:coord_dims]) - for target in instruction.targets_copy(): - coordinates[int(target.value)] = coord - return coordinates - - def _extract_mpp_layers(circuit: stim.Circuit) -> list[list[_MppProductRecord]]: layers: list[list[_MppProductRecord]] = [] current_layer: list[_MppProductRecord] | None = None @@ -250,7 +120,7 @@ def _extract_mpp_layers(circuit: stim.Circuit) -> list[list[_MppProductRecord]]: if instruction.name == "MPP": if current_layer is None: current_layer = [] - products = _mpp_targets_to_products(instruction.targets_copy()) + products = mpp_targets_to_products(instruction.targets_copy()) if len(products) != instruction.num_measurements: msg = "Stim MPP instruction measurement count does not match its parsed product count." raise ValueError(msg) @@ -281,14 +151,8 @@ def _select_mpp_products( return list(layers[mpp_layer]) -def _extract_selected_mpp_annotations( - circuit: stim.Circuit, - *, - record_to_row: Mapping[int, int], -) -> _StimMppAnnotations: - detector_rows: list[frozenset[int]] = [] +def _extract_stim_record_annotations(circuit: stim.Circuit) -> _StimRecordAnnotations: detector_record_indices: list[frozenset[int]] = [] - logical_observable_rows: dict[int, set[int]] = {} logical_observable_record_indices: dict[int, set[int]] = {} measurement_count = 0 @@ -298,138 +162,27 @@ def _extract_selected_mpp_annotations( raise TypeError(msg) if instruction.name == "DETECTOR": - rows, record_indices = _record_targets_to_selected_mpp_rows( + record_indices = record_targets_to_absolute_indices( instruction.targets_copy(), measurement_count=measurement_count, - record_to_row=record_to_row, instruction_name=instruction.name, ) - if rows is not None: - detector_rows.append(frozenset(rows)) - detector_record_indices.append(record_indices) + detector_record_indices.append(record_indices) elif instruction.name == "OBSERVABLE_INCLUDE": - logical_idx = _observable_index(instruction) - rows, record_indices = _record_targets_to_selected_mpp_rows( + logical_idx = observable_index(instruction) + record_indices = record_targets_to_absolute_indices( instruction.targets_copy(), measurement_count=measurement_count, - record_to_row=record_to_row, instruction_name=f"OBSERVABLE_INCLUDE({logical_idx})", ) - if rows is not None: - logical_observable_rows.setdefault(logical_idx, set()).symmetric_difference_update(rows) - logical_observable_record_indices.setdefault(logical_idx, set()).symmetric_difference_update( - record_indices - ) + logical_observable_record_indices.setdefault(logical_idx, set()).symmetric_difference_update(record_indices) measurement_count += instruction.num_measurements - return _StimMppAnnotations( - detector_rows=tuple(detector_rows), - logical_observable_rows={ - logical_idx: frozenset(rows) for logical_idx, rows in sorted(logical_observable_rows.items()) - }, + return _StimRecordAnnotations( detector_record_indices=tuple(detector_record_indices), logical_observable_record_indices={ logical_idx: frozenset(records) for logical_idx, records in sorted(logical_observable_record_indices.items()) }, ) - - -def _record_targets_to_selected_mpp_rows( - targets: Sequence[stim.GateTarget], - *, - measurement_count: int, - record_to_row: Mapping[int, int], - instruction_name: str, -) -> tuple[set[int] | None, frozenset[int]]: - rows: set[int] = set() - record_indices: set[int] = set() - saw_selected_record = False - - for target in targets: - if not target.is_measurement_record_target: - msg = f"{instruction_name} contains unsupported target {target!r}; only rec targets are supported." - raise ValueError(msg) - record_index = measurement_count + int(target.value) - if record_index in record_indices: - record_indices.remove(record_index) - else: - record_indices.add(record_index) - row = record_to_row.get(record_index) - if row is None: - continue - saw_selected_record = True - if row in rows: - rows.remove(row) - else: - rows.add(row) - - if not saw_selected_record: - return None, frozenset(record_indices) - return rows, frozenset(record_indices) - - -def _observable_index(instruction: stim.CircuitInstruction) -> int: - args = instruction.gate_args_copy() - if len(args) != 1 or not args[0].is_integer(): - msg = "OBSERVABLE_INCLUDE must have one integer observable index." - raise ValueError(msg) - return int(args[0]) - - -def _mpp_targets_to_products(targets: Sequence[stim.GateTarget]) -> list[PauliSupport]: - products: list[PauliSupport] = [] - current: list[tuple[int, str]] = [] - seen_in_current: set[int] = set() - expect_pauli = True - - for target in targets: - if target.is_combiner: - if expect_pauli: - msg = "Invalid MPP target list: unexpected combiner." - raise ValueError(msg) - expect_pauli = True - continue - - if target.is_inverted_result_target: - msg = "Signed MPP products are not supported; inverted Pauli targets cannot be imported." - raise ValueError(msg) - pauli = _target_pauli(target) - if current and not expect_pauli: - products.append(tuple(current)) - current = [] - seen_in_current = set() - - qid = int(target.value) - if qid in seen_in_current: - msg = f"Invalid MPP product: qubit {qid} appears more than once." - raise ValueError(msg) - current.append((qid, pauli)) - seen_in_current.add(qid) - expect_pauli = False - - if expect_pauli: - msg = "Invalid MPP target list: trailing combiner or empty product." - raise ValueError(msg) - products.append(tuple(current)) - return products - - -def _target_pauli(target: stim.GateTarget) -> str: - if target.is_x_target: - return "X" - if target.is_y_target: - return "Y" - if target.is_z_target: - return "Z" - msg = f"Unsupported MPP target: {target!r}." - raise ValueError(msg) - - -def _map_rows_to_nodes(rows: frozenset[int], ancilla_nodes: Mapping[int, int], label: str) -> set[int]: - missing_rows = sorted(row for row in rows if row not in ancilla_nodes) - if missing_rows: - msg = f"Cannot map {label}; ancilla node map is missing stabilizer row(s): {missing_rows}." - raise ValueError(msg) - return {ancilla_nodes[row] for row in rows} diff --git a/graphqomb/simulator.py b/graphqomb/simulator.py index 559bc2ec..1d523dbb 100644 --- a/graphqomb/simulator.py +++ b/graphqomb/simulator.py @@ -273,7 +273,11 @@ def simulate(self, rng: np.random.Generator | None = None) -> None: if node in self.__pattern.output_node_indices: self._apply_output_pauli_frame(node) - # Create a mapping from current node indices to output node indices - permutation = [self.__pattern.output_node_indices[node] for node in self.node_indices] + # Measured outputs can leave sparse qindices among the remaining quantum + # outputs. Reorder by each qindex's relative position, not by the qindex + # value itself. + output_qindices = [self.__pattern.output_node_indices[node] for node in self.node_indices] + qindex_rank = {qindex: rank for rank, qindex in enumerate(sorted(output_qindices))} + permutation = [qindex_rank[qindex] for qindex in output_qindices] self.state.reorder(permutation) diff --git a/graphqomb/statevec.py b/graphqomb/statevec.py index 1af2333d..ffa31bd1 100644 --- a/graphqomb/statevec.py +++ b/graphqomb/statevec.py @@ -55,8 +55,16 @@ def __init__(self, state: ArrayLike | None = None, *, copy: bool | None = None) # Internal qubit ordering: maps external qubit index to internal index self.__qindex_mng = QubitIndexManager(num_qubits) - def __array__(self, dtype: DTypeLike | None = None, copy: bool | None = None) -> NDArray[np.complex128]: - return np.asarray(self.state(), dtype=dtype, copy=copy) + def __array__(self, dtype: DTypeLike | None = None, copy: bool | None = None) -> NDArray[np.generic]: + state = self.state() + if dtype is not None: + target_dtype = np.dtype(dtype) + if target_dtype.kind in "bifu": + if np.any(state.imag != 0): + msg = "Cannot convert a state vector with nonzero imaginary amplitudes to a real dtype." + raise ValueError(msg) + return np.asarray(state.real, dtype=target_dtype, copy=copy) + return np.asarray(state, dtype=dtype, copy=copy) @property @typing_extensions.override diff --git a/graphqomb/stim_compiler.py b/graphqomb/stim_compiler.py index db62cd84..9b93b767 100644 --- a/graphqomb/stim_compiler.py +++ b/graphqomb/stim_compiler.py @@ -11,7 +11,7 @@ from typing import TYPE_CHECKING from graphqomb.command import TICK, E, M, N -from graphqomb.common import Axis, MeasBasis, determine_pauli_axis +from graphqomb.common import Axis, AxisMeasBasis, MeasBasis, Sign, determine_pauli_axis from graphqomb.noise_model import ( Coordinate, EntangleEvent, @@ -156,10 +156,12 @@ def _handle_measure(self, node: int, meas_basis: MeasBasis) -> None: # Emit measurement with optional flip probability meas_instr = {Axis.X: "MX", Axis.Y: "MY", Axis.Z: "MZ"}[axis] + is_inverted = isinstance(meas_basis, AxisMeasBasis) and meas_basis.sign is Sign.MINUS + meas_target = f"!{node}" if is_inverted else str(node) if meas_flip_p > 0.0: - self._stim_io.write(f"{meas_instr}({meas_flip_p}) {node}\n") + self._stim_io.write(f"{meas_instr}({meas_flip_p}) {meas_target}\n") else: - self._stim_io.write(f"{meas_instr} {node}\n") + self._stim_io.write(f"{meas_instr} {meas_target}\n") self._meas_order[node] = self._rec_index self._rec_index += 1 diff --git a/graphqomb/stim_importer.py b/graphqomb/stim_importer.py new file mode 100644 index 00000000..1b64966d --- /dev/null +++ b/graphqomb/stim_importer.py @@ -0,0 +1,859 @@ +"""Import supported Stim circuits into GraphQOMB patterns.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from graphlib import CycleError, TopologicalSorter +from itertools import combinations, pairwise +from pathlib import Path +from typing import TYPE_CHECKING + +import stim + +from graphqomb.circuit import Circuit, CircuitScheduleStrategy, circuit2graph +from graphqomb.common import Axis, AxisMeasBasis, Sign +from graphqomb.feedforward import dag_from_flow +from graphqomb.gates import CNOT, CZ, SWAP, Gate, H, Rz, S, X, Y, Z +from graphqomb.graphstate import GraphState, compose +from graphqomb.qec._stim import ( + PauliSupport, + StimMppExtraction, + extract_qubit_coordinates, + mpp_targets_to_products, + observable_index, + pauli_products_commute, + plain_qubit_target, + record_targets_to_absolute_indices, + stim_mpp_extraction_from_records, +) +from graphqomb.qec.qeccode import StabilizerGraphStateBuildResult, YFoliation, build_graph_state +from graphqomb.qompiler import qompile + +if TYPE_CHECKING: + from collections.abc import Callable, Mapping, Sequence + + from graphqomb.pattern import Pattern + + +_UNITARY_GATES = frozenset({"H", "S", "SQRT_Z", "S_DAG", "SQRT_Z_DAG", "X", "Y", "Z", "CX", "CNOT", "CZ", "SWAP"}) +_SINGLE_PAULI_MEASUREMENT_AXES = {"M": Axis.Z, "MX": Axis.X, "MY": Axis.Y} +_PAIR_PAULI_MEASUREMENT_AXES = {"MXX": "X", "MYY": "Y", "MZZ": "Z"} +_PAULI_PRODUCT_MEASUREMENT_GATES = frozenset({"MPP", *_PAIR_PAULI_MEASUREMENT_AXES}) +_SINGLE_QUBIT_GATE_FACTORIES: dict[str, Callable[[int], Gate]] = { + "H": H, + "S": S, + "SQRT_Z": S, + "S_DAG": lambda qubit: Rz(qubit, -math.pi / 2), + "SQRT_Z_DAG": lambda qubit: Rz(qubit, -math.pi / 2), + "X": X, + "Y": Y, + "Z": Z, +} +_TWO_QUBIT_GATE_FACTORIES: dict[str, Callable[[tuple[int, int]], Gate]] = { + "CX": CNOT, + "CNOT": CNOT, + "CZ": CZ, + "SWAP": SWAP, +} + + +@dataclass(frozen=True) +class StimImportResult: + """Result of importing a supported Stim circuit.""" + + pattern: Pattern + stim_to_qubit: dict[int, int] + qubit_to_stim: dict[int, int] + mpp_extractions: tuple[StimMppExtraction, ...] + + +@dataclass(frozen=True) +class _Fragment: + graph: GraphState + xflow: dict[int, set[int]] + record_nodes: dict[int, int] + mpp_extractions: tuple[StimMppExtraction, ...] = () + + +@dataclass(frozen=True) +class _ImportContext: + stim_to_qubit: Mapping[int, int] + coordinate_by_stim_id: Mapping[int, tuple[float, ...]] + detector_record_indices: Sequence[frozenset[int]] + logical_observable_record_indices: Mapping[int, frozenset[int]] + schedule_strategy: CircuitScheduleStrategy + y_foliation: YFoliation + + +@dataclass(frozen=True) +class _IdealizedCircuit: + circuit: stim.Circuit + zero_record_indices: frozenset[int] + + +@dataclass(frozen=True) +class _AnalyzedInstruction: + instruction: stim.CircuitInstruction + record_indices: tuple[int, ...] + + +@dataclass(frozen=True) +class _CircuitAnalysis: + blocks: tuple[tuple[_AnalyzedInstruction, ...], ...] + detector_record_indices: tuple[frozenset[int], ...] + logical_observable_record_indices: dict[int, frozenset[int]] + measurement_count: int + + +def stim_file_to_pattern( + path: str | Path, + *, + coord_dims: int = 2, + schedule_strategy: CircuitScheduleStrategy = CircuitScheduleStrategy.PARALLEL, + y_foliation: YFoliation = YFoliation.TYPE_I, +) -> StimImportResult: + """Import a supported Stim file into a GraphQOMB pattern. + + Returns + ------- + `StimImportResult` + Imported pattern, qubit mapping, and MPP extraction metadata. + """ + return stim_text_to_pattern( + Path(path).read_text(encoding="utf-8"), + coord_dims=coord_dims, + schedule_strategy=schedule_strategy, + y_foliation=y_foliation, + ) + + +def stim_text_to_pattern( + text: str, + *, + coord_dims: int = 2, + schedule_strategy: CircuitScheduleStrategy = CircuitScheduleStrategy.PARALLEL, + y_foliation: YFoliation = YFoliation.TYPE_I, +) -> StimImportResult: + """Import supported Stim text into a GraphQOMB pattern. + + Returns + ------- + `StimImportResult` + Imported pattern, qubit mapping, and MPP extraction metadata. + """ + return stim_circuit_to_pattern( + stim.Circuit(text), + coord_dims=coord_dims, + schedule_strategy=schedule_strategy, + y_foliation=y_foliation, + ) + + +def stim_circuit_to_pattern( + circuit: stim.Circuit, + *, + coord_dims: int = 2, + schedule_strategy: CircuitScheduleStrategy = CircuitScheduleStrategy.PARALLEL, + y_foliation: YFoliation = YFoliation.TYPE_I, +) -> StimImportResult: + """Import a supported Stim circuit into a GraphQOMB pattern. + + The importer supports Clifford unitary blocks and Pauli measurement blocks. + Stim noise instructions and measurement-error probabilities are omitted + because circuit-level noise is outside the GraphQOMB import model. Pauli + measurement blocks must be separated from unitary blocks by TICK. A direct + single-qubit measurement terminates that qubit's lifetime; other qubits may + continue, but the measured qubit cannot be used by a later operation. + + Returns + ------- + `StimImportResult` + Imported pattern, qubit mapping, and MPP extraction metadata. + + Raises + ------ + ValueError + If the circuit uses unsupported instructions or invalid coordinates. + """ + if coord_dims not in {2, 3}: + msg = "coord_dims must be 2 or 3." + raise ValueError(msg) + + flat_circuit = circuit.flattened() + idealized = _idealize_circuit(flat_circuit) + analysis = _analyze_circuit(idealized.circuit) + if analysis.measurement_count == 0 and ( + analysis.detector_record_indices or analysis.logical_observable_record_indices + ): + msg = "DETECTOR and OBSERVABLE_INCLUDE require at least one imported measurement instruction." + raise ValueError(msg) + coordinate_by_stim_id = extract_qubit_coordinates(idealized.circuit, coord_dims=coord_dims) + stim_to_qubit = _stim_to_qubit_map(idealized.circuit) + qubit_to_stim = {qubit: stim_id for stim_id, qubit in stim_to_qubit.items()} + context = _ImportContext( + stim_to_qubit=stim_to_qubit, + coordinate_by_stim_id=coordinate_by_stim_id, + detector_record_indices=analysis.detector_record_indices, + logical_observable_record_indices=analysis.logical_observable_record_indices, + schedule_strategy=schedule_strategy, + y_foliation=y_foliation, + ) + + fragments = _fragments_from_blocks(analysis.blocks, context=context) + fragment = _compose_fragments(fragments) + parity_check_groups, logical_observables = _measurement_annotations_from_analysis( + analysis, + record_nodes=fragment.record_nodes, + zero_record_indices=idealized.zero_record_indices, + ) + pattern = qompile( + fragment.graph, + fragment.xflow, + None, + parity_check_group=parity_check_groups, + logical_observables=logical_observables, + ) + return StimImportResult( + pattern=pattern, + stim_to_qubit=stim_to_qubit, + qubit_to_stim=qubit_to_stim, + mpp_extractions=fragment.mpp_extractions, + ) + + +def _idealize_circuit(circuit: stim.Circuit) -> _IdealizedCircuit: + """Remove circuit-level noise while preserving ideal measurement records. + + Returns + ------- + `_IdealizedCircuit` + Ideal circuit and absolute indices of zero-valued records. + + Raises + ------ + TypeError + If a flattened circuit unexpectedly contains a repeat block. + ValueError + If a constant true record cannot be represented. + """ + result = stim.Circuit() + zero_record_indices: set[int] = set() + measurement_offset = 0 + + for instruction in circuit: + if not isinstance(instruction, stim.CircuitInstruction): + msg = "Flattened Stim circuit unexpectedly contains a repeat block." + raise TypeError(msg) + + gate_data = stim.gate_data(instruction.name) + if instruction.name in _SINGLE_PAULI_MEASUREMENT_AXES: + result.append(instruction.name, instruction.targets_copy()) + elif instruction.name in _PAULI_PRODUCT_MEASUREMENT_GATES: + _append_ideal_pauli_measurements(result, instruction) + elif instruction.name == "MPAD": + targets = instruction.targets_copy() + if any(int(target.value) != 0 for target in targets): + msg = "MPAD 1 records are not supported because detector parity offsets are not represented." + raise ValueError(msg) + result.append("MPAD", targets) + zero_record_indices.update(range(measurement_offset, measurement_offset + instruction.num_measurements)) + elif gate_data.is_noisy_gate and not gate_data.is_reset: + if gate_data.produces_measurements: + result.append("MPAD", [0] * instruction.num_measurements) + zero_record_indices.update(range(measurement_offset, measurement_offset + instruction.num_measurements)) + else: + result.append(instruction) + + measurement_offset += instruction.num_measurements + + return _IdealizedCircuit(result, frozenset(zero_record_indices)) + + +def _analyze_circuit(circuit: stim.Circuit) -> _CircuitAnalysis: + """Index measurement records and split a flattened circuit at TICKs. + + Returns + ------- + `_CircuitAnalysis` + TICK-separated instructions and whole-circuit record annotations. + + Raises + ------ + TypeError + If a flattened circuit unexpectedly contains a repeat block. + """ + blocks: list[tuple[_AnalyzedInstruction, ...]] = [] + current_block: list[_AnalyzedInstruction] = [] + detector_record_indices: list[frozenset[int]] = [] + logical_record_indices: dict[int, set[int]] = {} + measurement_count = 0 + + for instruction in circuit: + if not isinstance(instruction, stim.CircuitInstruction): + msg = "Flattened Stim circuit unexpectedly contains a repeat block." + raise TypeError(msg) + if instruction.name == "TICK": + blocks.append(tuple(current_block)) + current_block = [] + continue + if instruction.name == "QUBIT_COORDS": + continue + if instruction.name == "DETECTOR": + detector_record_indices.append( + record_targets_to_absolute_indices( + instruction.targets_copy(), + measurement_count=measurement_count, + instruction_name=instruction.name, + ) + ) + elif instruction.name == "OBSERVABLE_INCLUDE": + logical_idx = observable_index(instruction) + logical_record_indices.setdefault(logical_idx, set()).symmetric_difference_update( + record_targets_to_absolute_indices( + instruction.targets_copy(), + measurement_count=measurement_count, + instruction_name=f"OBSERVABLE_INCLUDE({logical_idx})", + ) + ) + elif instruction.name != "MPAD": + record_indices = tuple(range(measurement_count, measurement_count + instruction.num_measurements)) + current_block.append(_AnalyzedInstruction(instruction, record_indices)) + + measurement_count += instruction.num_measurements + + blocks.append(tuple(current_block)) + return _CircuitAnalysis( + blocks=tuple(blocks), + detector_record_indices=tuple(detector_record_indices), + logical_observable_record_indices={ + logical_idx: frozenset(records) for logical_idx, records in sorted(logical_record_indices.items()) + }, + measurement_count=measurement_count, + ) + + +def _append_ideal_pauli_measurements( + circuit: stim.Circuit, + instruction: stim.CircuitInstruction, +) -> None: + """Append ideal MPP equivalents of Stim Pauli measurements. + + Raises + ------ + ValueError + If an instruction has an invalid target group. + """ + if instruction.name == "MPP": + circuit.append("MPP", instruction.targets_copy()) + return + + axis = _PAIR_PAULI_MEASUREMENT_AXES[instruction.name] + expected_group_size = 2 + + target_factory = {"X": stim.target_x, "Y": stim.target_y, "Z": stim.target_z}[axis] + for group in instruction.target_groups(): + if len(group) != expected_group_size: + msg = f"{instruction.name} target group must contain {expected_group_size} qubit(s)." + raise ValueError(msg) + product_targets: list[stim.GateTarget] = [] + for index, target in enumerate(group): + if index: + product_targets.append(stim.target_combiner()) + product_targets.append( + target_factory( + plain_qubit_target(target, instruction.name), + invert=target.is_inverted_result_target, + ) + ) + circuit.append("MPP", product_targets) + + +def _fragments_from_blocks( + blocks: Sequence[Sequence[_AnalyzedInstruction]], + *, + context: _ImportContext, +) -> list[_Fragment]: + _validate_blocks(blocks) + _validate_single_measurement_lifetimes(blocks) + + fragments = [_identity_fragment(context)] + mpp_layer_index = 0 + for block in blocks: + unitary_instructions = tuple( + analyzed.instruction for analyzed in block if analyzed.instruction.name in _UNITARY_GATES + ) + if unitary_instructions: + fragments.append(_unitary_fragment(unitary_instructions, context=context)) + else: + measurement_fragments = _measurement_fragments_from_block( + block, + mpp_layer_index=mpp_layer_index, + context=context, + ) + fragments.extend(measurement_fragments) + if any(analyzed.instruction.name == "MPP" for analyzed in block): + mpp_layer_index += sum( + len(analyzed.record_indices) for analyzed in block if analyzed.instruction.name == "MPP" + ) + + return fragments + + +def _validate_blocks(blocks: Sequence[Sequence[_AnalyzedInstruction]]) -> None: + """Validate supported instructions and required TICK separation. + + Raises + ------ + ValueError + If an instruction is unsupported or a block mixes unitary gates with + Pauli measurements. + """ + for block in blocks: + unsupported = [ + analyzed.instruction.name + for analyzed in block + if ( + analyzed.instruction.name not in _UNITARY_GATES + and analyzed.instruction.name not in _SINGLE_PAULI_MEASUREMENT_AXES + and analyzed.instruction.name != "MPP" + ) + ] + if unsupported: + msg = f"Unsupported Stim instruction(s): {', '.join(sorted(set(unsupported)))}." + raise ValueError(msg) + + has_unitary = any(analyzed.instruction.name in _UNITARY_GATES for analyzed in block) + has_pauli_measurement = any( + analyzed.instruction.name == "MPP" or analyzed.instruction.name in _SINGLE_PAULI_MEASUREMENT_AXES + for analyzed in block + ) + if has_unitary and has_pauli_measurement: + msg = "Pauli measurement instructions must be separated from unitary gate instructions by TICK." + raise ValueError(msg) + + +def _validate_single_measurement_lifetimes( + blocks: Sequence[Sequence[_AnalyzedInstruction]], +) -> None: + """Reject quantum operations after a directly measured qubit terminates. + + Raises + ------ + ValueError + If a quantum operation reuses a directly measured qubit. + """ + measured_qubits: set[int] = set() + for block in blocks: + for analyzed in block: + instruction = analyzed.instruction + if ( + instruction.name not in _UNITARY_GATES + and instruction.name != "MPP" + and instruction.name not in _SINGLE_PAULI_MEASUREMENT_AXES + ): + continue + + instruction_qubits = { + int(target.qubit_value) for target in instruction.targets_copy() if target.qubit_value is not None + } + reused_qubits = measured_qubits & instruction_qubits + if reused_qubits: + msg = ( + f"Stim qubit(s) {sorted(reused_qubits)} are used after a single-qubit measurement; " + "single-qubit measurements terminate those qubit lifetimes." + ) + raise ValueError(msg) + if instruction.name in _SINGLE_PAULI_MEASUREMENT_AXES: + measured_qubits.update(instruction_qubits) + + +def _measurement_fragments_from_block( + block: Sequence[_AnalyzedInstruction], + *, + mpp_layer_index: int, + context: _ImportContext, +) -> list[_Fragment]: + """Build direct and product measurement fragments in record order. + + Returns + ------- + `list`[`_Fragment`] + Measurement fragments in source order. + """ + fragments: list[_Fragment] = [] + mpp_items = tuple(analyzed for analyzed in block if analyzed.instruction.name == "MPP") + mpp_added = False + for analyzed in block: + instruction = analyzed.instruction + if instruction.name == "MPP" and not mpp_added: + fragments.append( + _mpp_fragment( + mpp_items, + mpp_layer_index=mpp_layer_index, + context=context, + ) + ) + mpp_added = True + elif instruction.name in _SINGLE_PAULI_MEASUREMENT_AXES: + fragments.append( + _single_measurement_fragment( + instruction, + record_indices=analyzed.record_indices, + context=context, + ) + ) + return fragments + + +def _single_measurement_fragment( + instruction: stim.CircuitInstruction, + *, + record_indices: Sequence[int], + context: _ImportContext, +) -> _Fragment: + """Build a fragment by assigning a basis directly to each measured node. + + Returns + ------- + `_Fragment` + Direct-measurement graph fragment and record-to-node mapping. + + Raises + ------ + ValueError + If target counts differ or a qubit is repeated in the instruction. + """ + targets = instruction.targets_copy() + if len(targets) != len(record_indices): + msg = f"{instruction.name} target count does not match its measurement-record count." + raise ValueError(msg) + + graph = GraphState() + record_nodes: dict[int, int] = {} + seen_qubits: set[int] = set() + axis = _SINGLE_PAULI_MEASUREMENT_AXES[instruction.name] + for target, record_index in zip(targets, record_indices, strict=True): + stim_id = plain_qubit_target(target, instruction.name) + if stim_id in seen_qubits: + msg = f"{instruction.name} measures qubit {stim_id} more than once in one instruction." + raise ValueError(msg) + seen_qubits.add(stim_id) + + node = graph.add_node(coordinate=context.coordinate_by_stim_id.get(stim_id)) + qubit_index = context.stim_to_qubit[stim_id] + graph.register_input(node, qubit_index) + graph.register_output(node, qubit_index) + sign = Sign.MINUS if target.is_inverted_result_target else Sign.PLUS + graph.assign_meas_basis(node, AxisMeasBasis(axis, sign)) + record_nodes[record_index] = node + + return _Fragment( + graph=graph, + xflow={}, + record_nodes=record_nodes, + ) + + +def _identity_fragment(context: _ImportContext) -> _Fragment: + graph = GraphState() + for stim_id, qubit_index in sorted(context.stim_to_qubit.items()): + node = graph.add_node(coordinate=context.coordinate_by_stim_id.get(stim_id)) + graph.register_input(node, qubit_index) + graph.register_output(node, qubit_index) + return _Fragment(graph=graph, xflow={}, record_nodes={}) + + +def _unitary_fragment( + block: Sequence[stim.CircuitInstruction], + *, + context: _ImportContext, +) -> _Fragment: + active_stim_ids = sorted( + {plain_qubit_target(target, instruction.name) for instruction in block for target in instruction.targets_copy()} + ) + stim_to_local = {stim_id: local_index for local_index, stim_id in enumerate(active_stim_ids)} + local_to_global = {local_index: context.stim_to_qubit[stim_id] for stim_id, local_index in stim_to_local.items()} + circuit = Circuit(len(active_stim_ids)) + for instruction in block: + _append_unitary_instruction(circuit, instruction, stim_to_local) + + local_graph, local_xflow, _scheduler = circuit2graph(circuit, schedule_strategy=context.schedule_strategy) + graph, node_map = _copy_graph_with_qindices(local_graph, local_to_global) + _apply_stim_coordinates( + graph, + stim_to_qubit=context.stim_to_qubit, + coordinate_by_stim_id=context.coordinate_by_stim_id, + ) + return _Fragment( + graph=graph, + xflow=_remap_flow(local_xflow, node_map), + record_nodes={}, + ) + + +def _copy_graph_with_qindices( + graph: GraphState, + local_to_global: Mapping[int, int], +) -> tuple[GraphState, dict[int, int]]: + """Copy a graph while replacing its input and output qindices. + + Returns + ------- + tuple[GraphState, dict[int, int]] + Copied graph and source-to-copy node map. + """ + copied = GraphState() + node_map = {node: copied.add_node() for node in sorted(graph.nodes)} + for node1, node2 in graph.edges: + copied.add_edge(node_map[node1], node_map[node2]) + for node, q_index in graph.input_node_indices.items(): + copied.register_input(node_map[node], local_to_global[q_index]) + for node, q_index in graph.output_node_indices.items(): + copied.register_output(node_map[node], local_to_global[q_index]) + for node, meas_basis in graph.meas_bases.items(): + copied.assign_meas_basis(node_map[node], meas_basis) + for node, local_clifford in graph.local_cliffords.items(): + copied.apply_local_clifford(node_map[node], local_clifford) + for node, coordinate in graph.coordinates.items(): + copied.set_coordinate(node_map[node], coordinate) + return copied, node_map + + +def _append_unitary_instruction( + circuit: Circuit, + instruction: stim.CircuitInstruction, + stim_to_qubit: Mapping[int, int], +) -> None: + for group in instruction.target_groups(): + qubits = [plain_qubit_target(target, instruction.name) for target in group] + mapped = [stim_to_qubit[qubit] for qubit in qubits] + single_factory = _SINGLE_QUBIT_GATE_FACTORIES.get(instruction.name) + two_factory = _TWO_QUBIT_GATE_FACTORIES.get(instruction.name) + if single_factory is not None: + circuit.apply_macro_gate(single_factory(mapped[0])) + elif two_factory is not None: + circuit.apply_macro_gate(two_factory((mapped[0], mapped[1]))) + else: + msg = f"Unsupported unitary Stim instruction: {instruction.name}." + raise ValueError(msg) + + +def _mpp_fragment( + block: Sequence[_AnalyzedInstruction], + *, + mpp_layer_index: int, + context: _ImportContext, +) -> _Fragment: + supports = tuple( + support for analyzed in block for support in mpp_targets_to_products(analyzed.instruction.targets_copy()) + ) + _validate_commuting_mpp_supports(supports) + record_indices = tuple(record_index for analyzed in block for record_index in analyzed.record_indices) + extraction = stim_mpp_extraction_from_records( + supports, + record_indices, + coordinate_by_stim_id=context.coordinate_by_stim_id, + detector_record_indices=context.detector_record_indices, + logical_observable_record_indices=context.logical_observable_record_indices, + ) + z_base = 2 * mpp_layer_index + fragment = _mpp_graph_fragment( + extraction, + record_indices=record_indices, + z_base=z_base, + context=context, + ) + if _has_causal_flow(fragment): + return _with_mpp_extraction(fragment, extraction) + + serialized_fragments = [ + _mpp_graph_fragment( + stim_mpp_extraction_from_records( + (support,), + (record_index,), + coordinate_by_stim_id=context.coordinate_by_stim_id, + detector_record_indices=context.detector_record_indices, + logical_observable_record_indices=context.logical_observable_record_indices, + ), + record_indices=(record_index,), + z_base=z_base + 2 * row, + context=context, + ) + for row, (support, record_index) in enumerate(zip(supports, record_indices, strict=True)) + ] + return _with_mpp_extraction(_compose_fragments(serialized_fragments), extraction) + + +def _validate_commuting_mpp_supports(supports: Sequence[PauliSupport]) -> None: + for (left_index, left), (right_index, right) in combinations(enumerate(supports), 2): + if not pauli_products_commute(left, right): + msg = ( + f"MPP products within one TICK block must commute; products {left_index} and {right_index} anticommute." + ) + raise ValueError(msg) + + +def _mpp_graph_fragment( + extraction: StimMppExtraction, + *, + record_indices: Sequence[int], + z_base: int, + context: _ImportContext, +) -> _Fragment: + qubit_indices = {column: context.stim_to_qubit[stim_id] for column, stim_id in extraction.column_to_stim.items()} + result = build_graph_state( + extraction.code, + z_base=z_base, + y_foliation=context.y_foliation, + data_as_io=True, + qubit_indices=qubit_indices, + ) + xflow = _mpp_flow(result) + if len(record_indices) != len(result.ancilla_nodes): + msg = "Imported MPP record count does not match the generated ancilla-node count." + raise ValueError(msg) + return _Fragment( + graph=result.graph, + xflow=xflow, + record_nodes={record_indices[row]: node for row, node in result.ancilla_nodes.items()}, + ) + + +def _has_causal_flow(fragment: _Fragment) -> bool: + try: + tuple(TopologicalSorter(dag_from_flow(fragment.graph, fragment.xflow)).static_order()) + except CycleError: + return False + return True + + +def _with_mpp_extraction(fragment: _Fragment, extraction: StimMppExtraction) -> _Fragment: + return _Fragment( + graph=fragment.graph, + xflow=fragment.xflow, + record_nodes=fragment.record_nodes, + mpp_extractions=(extraction,), + ) + + +def _mpp_flow( + result: StabilizerGraphStateBuildResult, +) -> dict[int, set[int]]: + xflow: dict[int, set[int]] = {} + measured_nodes_by_qubit: dict[int, list[int]] = {} + for qubit in sorted({key[0] for key in result.data_nodes}): + layer_nodes = [node for (data_qubit, _layer), node in sorted(result.data_nodes.items()) if data_qubit == qubit] + measured_nodes_by_qubit[qubit] = [node for node in layer_nodes if node in result.graph.meas_bases] + for current_node, next_node in pairwise(layer_nodes): + if current_node in result.graph.meas_bases: + xflow[current_node] = {next_node} + for ancilla_node in result.ancilla_nodes.values(): + correction_nodes = {ancilla_node} + for measured_nodes in measured_nodes_by_qubit.values(): + for earlier_node, later_node in pairwise(measured_nodes): + # Type I Y support touches both data-measurement layers. Including + # the later data stabilizer cancels the backward dependency in + # the automatically derived odd-neighborhood zflow. + if result.graph.has_edge(ancilla_node, earlier_node) and result.graph.has_edge( + ancilla_node, later_node + ): + correction_nodes.add(later_node) + xflow[ancilla_node] = correction_nodes + return xflow + + +def _compose_fragments(fragments: Sequence[_Fragment]) -> _Fragment: + current = fragments[0] + for fragment in fragments[1:]: + graph, node_map1, node_map2 = compose(current.graph, fragment.graph) + current = _Fragment( + graph=graph, + xflow=_remap_flow(current.xflow, node_map1) | _remap_flow(fragment.xflow, node_map2), + record_nodes=_remap_record_nodes(current.record_nodes, node_map1) + | _remap_record_nodes(fragment.record_nodes, node_map2), + mpp_extractions=(*current.mpp_extractions, *fragment.mpp_extractions), + ) + return current + + +def _measurement_annotations_from_analysis( + analysis: _CircuitAnalysis, + *, + record_nodes: Mapping[int, int], + zero_record_indices: frozenset[int], +) -> tuple[list[set[int]], dict[int, set[int]]]: + if not record_nodes and not zero_record_indices: + return [], {} + + parity_check_groups = [ + _record_indices_to_nodes(record_indices, record_nodes, zero_record_indices=zero_record_indices) + for record_indices in analysis.detector_record_indices + ] + logical_observables = { + logical_idx: _record_indices_to_nodes( + record_indices, + record_nodes, + zero_record_indices=zero_record_indices, + ) + for logical_idx, record_indices in analysis.logical_observable_record_indices.items() + } + return parity_check_groups, logical_observables + + +def _record_indices_to_nodes( + record_indices: frozenset[int], + record_nodes: Mapping[int, int], + *, + zero_record_indices: frozenset[int], +) -> set[int]: + missing_records = sorted( + record_index + for record_index in record_indices + if record_index not in record_nodes and record_index not in zero_record_indices + ) + if missing_records: + msg = f"Cannot map Stim measurement record(s) to imported Pauli-measurement nodes: {missing_records}." + raise ValueError(msg) + return {record_nodes[record_index] for record_index in record_indices if record_index in record_nodes} + + +def _remap_flow(flow: Mapping[int, set[int]], node_map: Mapping[int, int]) -> dict[int, set[int]]: + return {node_map[node]: _remap_node_set(targets, node_map) for node, targets in flow.items()} + + +def _remap_node_set(nodes: set[int], node_map: Mapping[int, int]) -> set[int]: + return {node_map[node] for node in nodes} + + +def _remap_record_nodes(record_nodes: Mapping[int, int], node_map: Mapping[int, int]) -> dict[int, int]: + return {record_index: node_map[node] for record_index, node in record_nodes.items()} + + +def _apply_stim_coordinates( + graph: GraphState, + *, + stim_to_qubit: Mapping[int, int], + coordinate_by_stim_id: Mapping[int, tuple[float, ...]], +) -> None: + qubit_to_stim = {qubit: stim_id for stim_id, qubit in stim_to_qubit.items()} + for node, q_index in graph.input_node_indices.items() | graph.output_node_indices.items(): + stim_id = qubit_to_stim[q_index] + coord = coordinate_by_stim_id.get(stim_id) + if coord is not None: + graph.set_coordinate(node, coord) + + +def _stim_to_qubit_map(circuit: stim.Circuit) -> dict[int, int]: + stim_ids: set[int] = set() + for instruction in circuit: + if not isinstance(instruction, stim.CircuitInstruction): + msg = "Flattened Stim circuit unexpectedly contains a repeat block." + raise TypeError(msg) + if instruction.name in {"TICK", "DETECTOR", "OBSERVABLE_INCLUDE", "MPAD"}: + continue + for target in instruction.targets_copy(): + qubit_value = target.qubit_value + if qubit_value is not None: + stim_ids.add(int(qubit_value)) + return {stim_id: qubit for qubit, stim_id in enumerate(sorted(stim_ids))} diff --git a/tests/test_qec.py b/tests/test_qec.py index 06b42905..1d2b6bb7 100644 --- a/tests/test_qec.py +++ b/tests/test_qec.py @@ -74,6 +74,72 @@ def test_build_graph_state_returns_index_to_node_maps() -> None: assert set(result.data_nodes.values()).isdisjoint(result.ancilla_nodes.values()) +def test_build_graph_state_registers_custom_data_io_indices() -> None: + code = StabilizerCode(_matrix([[1, 0, 0, 1]])) + + result = build_graph_state(code, data_as_io=True, qubit_indices={0: 10, 1: 12}) + + assert result.graph.input_node_indices == {result.data_nodes[0, 0]: 10, result.data_nodes[1, 0]: 12} + assert result.graph.output_node_indices == {result.data_nodes[0, 2]: 10, result.data_nodes[1, 2]: 12} + assert set(result.data_nodes) == {(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)} + for qubit in range(2): + _assert_axis_meas_basis(result.graph.meas_bases[result.data_nodes[qubit, 0]], Axis.X) + _assert_axis_meas_basis(result.graph.meas_bases[result.data_nodes[qubit, 1]], Axis.X) + assert result.data_nodes[qubit, 2] not in result.graph.meas_bases + assert result.graph.has_edge(result.data_nodes[qubit, 0], result.data_nodes[qubit, 1]) + assert result.graph.has_edge(result.data_nodes[qubit, 1], result.data_nodes[qubit, 2]) + + +@pytest.mark.parametrize( + ("qubit_indices", "message"), + [ + ({0: 10}, "map every stabilizer-code qubit"), + ({0: 10, 1: 10}, "values must be unique"), + ], +) +def test_build_graph_state_rejects_invalid_custom_data_io_indices( + qubit_indices: dict[int, int], + message: str, +) -> None: + code = StabilizerCode(_matrix([[1, 0, 0, 1]])) + + with pytest.raises(ValueError, match=message): + build_graph_state(code, data_as_io=True, qubit_indices=qubit_indices) + + +def test_build_graph_state_rejects_qubit_indices_without_data_io() -> None: + code = StabilizerCode(_matrix([[1, 0, 0, 1]])) + + with pytest.raises(ValueError, match="data_as_io=True"): + build_graph_state(code, qubit_indices={0: 10, 1: 12}) + + +def test_build_graph_state_registers_type_ii_chain_endpoints_as_io() -> None: + code = StabilizerCode(_matrix([[1, 1]])) + + result = build_graph_state(code, y_foliation=YFoliation.TYPE_II, data_as_io=True) + + assert result.graph.input_node_indices == {result.data_nodes[0, 0]: 0} + assert result.graph.output_node_indices == {result.data_nodes[0, 3]: 0} + assert set(result.data_nodes) == {(0, 0), (0, 1), (0, 2), (0, 3)} + for layer in range(3): + _assert_axis_meas_basis(result.graph.meas_bases[result.data_nodes[0, layer]], Axis.Y) + assert result.data_nodes[0, 3] not in result.graph.meas_bases + assert result.graph.has_edge(result.data_nodes[0, 2], result.data_nodes[0, 3]) + + +def test_build_graph_state_type_ii_keeps_separate_output_for_non_y_support() -> None: + code = StabilizerCode(_matrix([[1, 0]])) + + result = build_graph_state(code, y_foliation=YFoliation.TYPE_II, data_as_io=True) + + assert set(result.data_nodes) == {(0, 0), (0, 1), (0, 2)} + assert result.graph.output_node_indices == {result.data_nodes[0, 2]: 0} + _assert_axis_meas_basis(result.graph.meas_bases[result.data_nodes[0, 0]], Axis.X) + _assert_axis_meas_basis(result.graph.meas_bases[result.data_nodes[0, 1]], Axis.X) + assert result.data_nodes[0, 2] not in result.graph.meas_bases + + def test_build_graph_state_type_ii_routes_y_support_through_three_node_y_chain() -> None: code = StabilizerCode( _matrix( @@ -143,6 +209,17 @@ def test_build_graph_state_type_ii_aligns_three_node_chain_output_z_with_two_nod assert coords[result.ancilla_nodes[0]] == (10.0, 20.0, 5.5) +@pytest.mark.parametrize("y_foliation", [YFoliation.TYPE_I, YFoliation.TYPE_II]) +def test_build_graph_state_places_separate_io_output_at_end_of_unit(y_foliation: YFoliation) -> None: + code = StabilizerCode(_matrix([[1, 1]]), qubit_coords={0: (10.0, 20.0)}) + + result = build_graph_state(code, z_base=5, y_foliation=y_foliation, data_as_io=True) + output_node = next(iter(result.graph.output_node_indices)) + + assert result.graph.coordinates[output_node] == (10.0, 20.0, 7.0) + assert output_node not in result.graph.meas_bases + + def test_build_graph_state_lifts_coordinates_to_shifted_3d_layers() -> None: code = StabilizerCode( _matrix([[0, 1, 1, 0]]), diff --git a/tests/test_simulator.py b/tests/test_simulator.py index ec12f2fc..24313636 100644 --- a/tests/test_simulator.py +++ b/tests/test_simulator.py @@ -78,3 +78,52 @@ def test_pattern_simulator_measures_output_z_frame_in_x_basis() -> None: assert simulator.results == {node: True} assert simulator.output_results == {0: True} + + +def test_pattern_simulator_reorders_remaining_outputs_after_terminal_measurement() -> None: + """A measured output may leave sparse qindices for quantum outputs.""" + graph = GraphState() + measured_node = graph.add_node() + quantum_node = graph.add_node() + graph.register_input(measured_node, 0) + graph.register_input(quantum_node, 1) + graph.register_output(measured_node, 0) + graph.register_output(quantum_node, 1) + graph.assign_meas_basis(measured_node, AxisMeasBasis(Axis.Z, Sign.PLUS)) + pattern = Pattern( + input_node_indices=graph.input_node_indices, + output_node_indices=graph.output_node_indices, + commands=(M(measured_node, AxisMeasBasis(Axis.Z, Sign.PLUS)),), + pauli_frame=PauliFrame(graph, xflow={}, zflow={}), + ) + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + + simulator.simulate(rng=np.random.default_rng(3)) + + assert set(simulator.output_results) == {0} + np.testing.assert_allclose(simulator.state.state(), np.asarray([1, 1]) / np.sqrt(2)) + + +def test_pattern_simulator_reorders_outputs_by_qindex_rank() -> None: + """A non-self-inverse output permutation should preserve qindex ordering.""" + graph = GraphState() + for qindex in (2, 0, 1): + node = graph.add_node() + graph.register_input(node, qindex) + graph.register_output(node, qindex) + pattern = Pattern( + input_node_indices=graph.input_node_indices, + output_node_indices=graph.output_node_indices, + commands=(), + pauli_frame=PauliFrame(graph, xflow={}, zflow={}), + ) + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + initial_state = np.zeros((2, 2, 2), dtype=np.complex128) + initial_state[0, 1, 0] = 1 + simulator.state = StateVector(initial_state) + + simulator.simulate() + + expected_state = np.zeros((2, 2, 2), dtype=np.complex128) + expected_state[1, 0, 0] = 1 + np.testing.assert_allclose(simulator.state.state(), expected_state) diff --git a/tests/test_statevec.py b/tests/test_statevec.py index 5d821270..3e7e022c 100644 --- a/tests/test_statevec.py +++ b/tests/test_statevec.py @@ -515,11 +515,24 @@ def test_array_method_with_dtype() -> None: assert arr_complex64.dtype == np.complex64 assert np.allclose(arr_complex64.flatten(), state.astype(np.complex64)) - # Test with float64 dtype (should keep real part only) + # A state with real-valued amplitudes can be converted without warning arr_float = np.array(sv, dtype=np.float64) assert arr_float.dtype == np.float64 assert np.allclose(arr_float.flatten(), state.real.astype(np.float64)) + arr_int = np.array(sv, dtype=np.int64) + assert arr_int.dtype == np.int64 + assert np.array_equal(arr_int.flatten(), state.real.astype(np.int64)) + + +@pytest.mark.parametrize("dtype", [np.float64, np.int64, np.bool_]) +def test_array_method_rejects_noncomplex_dtype_for_complex_amplitudes(dtype: type[np.generic]) -> None: + """Test that non-complex conversion does not discard imaginary amplitudes.""" + sv = StateVector(np.array([1, 1j], dtype=np.complex128)) + + with pytest.raises(ValueError, match="nonzero imaginary amplitudes"): + np.array(sv, dtype=dtype) + def test_array_method_preserves_shape() -> None: """Test that __array__ method preserves the tensor structure for multi-qubit states.""" diff --git a/tests/test_stim_importer.py b/tests/test_stim_importer.py new file mode 100644 index 00000000..1f56758a --- /dev/null +++ b/tests/test_stim_importer.py @@ -0,0 +1,507 @@ +"""Tests for importing supported Stim circuits into GraphQOMB patterns.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +import pytest + +from graphqomb.command import M +from graphqomb.common import Axis, AxisMeasBasis, Sign +from graphqomb.graphstate import odd_neighbors +from graphqomb.qec.qeccode import YFoliation +from graphqomb.simulator import PatternSimulator, SimulatorBackend +from graphqomb.statevec import StateVector +from graphqomb.stim_compiler import stim_compile +from graphqomb.stim_importer import stim_circuit_to_pattern, stim_file_to_pattern, stim_text_to_pattern + +stim = pytest.importorskip("stim") + +if TYPE_CHECKING: + from pathlib import Path + + +def test_stim_text_to_pattern_imports_unitary_clifford_block() -> None: + result = stim_text_to_pattern( + """ + H 10 + CX 10 12 + S_DAG 12 + """ + ) + + assert result.stim_to_qubit == {10: 0, 12: 1} + assert result.qubit_to_stim == {0: 10, 1: 12} + assert result.mpp_extractions == () + assert set(result.pattern.input_node_indices.values()) == {0, 1} + assert set(result.pattern.output_node_indices.values()) == {0, 1} + + +def test_stim_text_to_pattern_preserves_unitary_semantics_across_ticks() -> None: + initial = np.asarray([1.0, 1.0], dtype=np.complex128) / np.sqrt(2) + expected = np.asarray([1 + 1j, 1 - 1j], dtype=np.complex128) / 2 + + for seed in range(8): + pattern = stim_text_to_pattern("S 0\nTICK\nH 0\n").pattern + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + simulator.state = StateVector(initial) + simulator.simulate(rng=np.random.default_rng(seed)) + + overlap = np.vdot(expected, simulator.state.state()) + assert np.isclose(abs(overlap), 1.0, atol=1e-9) + + +def test_stim_text_to_pattern_preserves_sparse_qubit_coordinates() -> None: + result = stim_text_to_pattern( + """ + QUBIT_COORDS(1, 2) 10 + QUBIT_COORDS(3) 99 + QUBIT_COORDS(4) 99 + H 99 + """ + ) + + assert result.stim_to_qubit == {10: 0, 99: 1} + assert result.pattern.input_coordinates + assert set(result.pattern.input_coordinates.values()) == {(1.0, 2.0), (3.0, 4.0)} + + +@pytest.mark.parametrize( + ("instruction", "initial", "expected"), + [ + ("H 0", [1, 0], [1 / np.sqrt(2), 1 / np.sqrt(2)]), + ("S 0", [1 / np.sqrt(2), 1 / np.sqrt(2)], [1 / np.sqrt(2), 1j / np.sqrt(2)]), + ("S_DAG 0", [1 / np.sqrt(2), 1 / np.sqrt(2)], [1 / np.sqrt(2), -1j / np.sqrt(2)]), + ("X 0", [1, 0], [0, 1]), + ("Y 0", [1, 0], [0, 1j]), + ("Z 0", [1 / np.sqrt(2), 1 / np.sqrt(2)], [1 / np.sqrt(2), -1 / np.sqrt(2)]), + ], +) +def test_stim_text_to_pattern_preserves_supported_single_qubit_gates( + instruction: str, + initial: list[complex], + expected: list[complex], +) -> None: + pattern = stim_text_to_pattern(instruction).pattern + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + simulator.state = StateVector(initial) + + simulator.simulate(rng=np.random.default_rng(3)) + + assert np.isclose(abs(np.vdot(expected, simulator.state.state())), 1.0, atol=1e-9) + + +@pytest.mark.parametrize( + ("instruction", "initial", "expected"), + [ + ("CX 0 1", [0, 0, 1, 0], [0, 0, 0, 1]), + ("CZ 0 1", [0.5, 0.5, 0.5, 0.5], [0.5, 0.5, 0.5, -0.5]), + ("SWAP 0 1", [0, 0, 1, 0], [0, 1, 0, 0]), + ], +) +def test_stim_text_to_pattern_preserves_supported_two_qubit_gates( + instruction: str, + initial: list[complex], + expected: list[complex], +) -> None: + pattern = stim_text_to_pattern(instruction).pattern + simulator = PatternSimulator(pattern, SimulatorBackend.StateVector) + simulator.state = StateVector(initial) + + simulator.simulate(rng=np.random.default_rng(3)) + + assert np.isclose(abs(np.vdot(expected, simulator.state.state())), 1.0, atol=1e-9) + + +def test_stim_text_to_pattern_imports_tick_separated_mpp_block() -> None: + result = stim_text_to_pattern( + """ + H 10 + TICK + MPP X10*Z12 + DETECTOR rec[-1] + OBSERVABLE_INCLUDE(3) rec[-1] + TICK + CZ 10 12 + """ + ) + + assert result.stim_to_qubit == {10: 0, 12: 1} + assert len(result.mpp_extractions) == 1 + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + assert set(result.pattern.pauli_frame.logical_observables) == {3} + assert set(result.pattern.output_node_indices.values()) == {0, 1} + + +def test_stim_text_to_pattern_combines_commuting_mpp_instructions_in_one_tick_block() -> None: + result = stim_text_to_pattern( + """ + MPP X0 + DETECTOR rec[-1] + MPP Z1 + DETECTOR rec[-1] + """ + ) + + assert len(result.mpp_extractions) == 1 + assert result.mpp_extractions[0].supports == (((0, "X"),), ((1, "Z"),)) + assert len(result.pattern.pauli_frame.parity_check_group) == 2 + + +def test_stim_text_to_pattern_rejects_anticommuting_mpp_in_one_tick_block() -> None: + with pytest.raises(ValueError, match="must commute"): + stim_text_to_pattern("MPP X0\nMPP Z0") + + +@pytest.mark.parametrize("y_foliation", [YFoliation.TYPE_I, YFoliation.TYPE_II]) +def test_stim_text_to_pattern_serializes_cyclic_commuting_mpp_flow(y_foliation: YFoliation) -> None: + result = stim_text_to_pattern( + """ + MPP X0*X1*X4*X5 + MPP Z0*Z1*Z2*Z3 + MPP Y0*X2*Z4*Z6 + MPP Z4*Z5 + MPP X1*X3 + MPP Z2*X6 + """, + y_foliation=y_foliation, + ) + + assert len(result.mpp_extractions) == 1 + assert len(result.mpp_extractions[0].supports) == 6 + assert set(result.pattern.input_node_indices.values()) == set(range(7)) + assert set(result.pattern.output_node_indices.values()) == set(range(7)) + + +def test_stim_text_to_pattern_uses_automatic_zflow_for_mpp_graph() -> None: + result = stim_text_to_pattern("MPP X0*Z1") + graph = result.pattern.pauli_frame.graphstate + + assert result.pattern.pauli_frame.xflow + for node, correction_nodes in result.pattern.pauli_frame.xflow.items(): + assert result.pattern.pauli_frame.zflow[node] == odd_neighbors(correction_nodes, graph) + + +def test_stim_text_to_pattern_appends_output_after_type_i_mpp_measurements() -> None: + result = stim_text_to_pattern("MPP X0") + graph = result.pattern.pauli_frame.graphstate + + assert graph.number_of_nodes() == 4 + assert len(graph.meas_bases) == 3 + assert len(graph.output_node_indices) == 1 + assert next(iter(graph.output_node_indices)) not in graph.meas_bases + + +def test_stim_text_to_pattern_appends_output_after_type_ii_y_measurements() -> None: + result = stim_text_to_pattern("MPP Y0", y_foliation=YFoliation.TYPE_II) + graph = result.pattern.pauli_frame.graphstate + y_measurements = [ + basis for basis in graph.meas_bases.values() if isinstance(basis, AxisMeasBasis) and basis.axis == Axis.Y + ] + + assert graph.number_of_nodes() == 5 + assert len(y_measurements) == 3 + assert len(graph.output_node_indices) == 1 + assert next(iter(graph.output_node_indices)) not in graph.meas_bases + + +def test_stim_import_entry_points_accept_type_ii_foliation(tmp_path: Path) -> None: + stim_path = tmp_path / "y_measurement.stim" + stim_path.write_text("MPP Y0", encoding="utf-8") + + circuit_result = stim_circuit_to_pattern(stim.Circuit("MPP Y0"), y_foliation=YFoliation.TYPE_II) + file_result = stim_file_to_pattern(stim_path, y_foliation=YFoliation.TYPE_II) + + for result in (circuit_result, file_result): + axes = [ + basis.axis + for basis in result.pattern.pauli_frame.graphstate.meas_bases.values() + if isinstance(basis, AxisMeasBasis) + ] + assert axes.count(Axis.Y) == 3 + + +@pytest.mark.parametrize( + ("instruction", "expected_axis"), + [ + ("M 10", Axis.Z), + ("MZ 10", Axis.Z), + ("MX 10", Axis.X), + ("MY 10", Axis.Y), + ], +) +def test_stim_text_to_pattern_imports_single_qubit_pauli_measurements( + instruction: str, + expected_axis: Axis, +) -> None: + result = stim_text_to_pattern(instruction) + + measurements = [command for command in result.pattern.commands if isinstance(command, M)] + + assert result.mpp_extractions == () + assert result.stim_to_qubit == {10: 0} + assert result.pattern.input_node_indices == {0: 0} + assert result.pattern.output_node_indices == {0: 0} + assert len(measurements) == 1 + assert measurements[0].node == 0 + assert isinstance(measurements[0].meas_basis, AxisMeasBasis) + assert measurements[0].meas_basis.axis == expected_axis + assert measurements[0].meas_basis.sign == Sign.PLUS + + +def test_stim_text_to_pattern_assigns_single_measurement_to_existing_wire_node() -> None: + result = stim_text_to_pattern("H 10\nTICK\nMX 10") + output_node = next(node for node, qubit in result.pattern.output_node_indices.items() if qubit == 0) + output_measurements = [ + command for command in result.pattern.commands if isinstance(command, M) and command.node == output_node + ] + + assert result.mpp_extractions == () + assert len(output_measurements) == 1 + assert isinstance(output_measurements[0].meas_basis, AxisMeasBasis) + assert output_measurements[0].meas_basis.axis == Axis.X + assert output_measurements[0].meas_basis.sign == Sign.PLUS + + +@pytest.mark.parametrize( + ("instruction", "expected_axis"), + [ + ("MXX 10 12", "X"), + ("MYY 10 12", "Y"), + ("MZZ 10 12", "Z"), + ], +) +def test_stim_text_to_pattern_imports_pair_pauli_measurements( + instruction: str, + expected_axis: str, +) -> None: + result = stim_text_to_pattern(instruction) + + assert result.mpp_extractions[0].supports == (((10, expected_axis), (12, expected_axis)),) + assert result.stim_to_qubit == {10: 0, 12: 1} + + +def test_stim_text_to_pattern_preserves_multiple_measurement_results_in_target_order() -> None: + result = stim_text_to_pattern( + """ + M 0 2 + MXX 1 3 4 5 + DETECTOR rec[-4] rec[-1] + """ + ) + + assert result.mpp_extractions[0].supports == ( + ((1, "X"), (3, "X")), + ((4, "X"), (5, "X")), + ) + direct_measurements = [command for command in result.pattern.commands if isinstance(command, M)] + assert ( + sum( + isinstance(command.meas_basis, AxisMeasBasis) and command.meas_basis.axis == Axis.Z + for command in direct_measurements + ) + == 2 + ) + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + assert len(result.pattern.pauli_frame.parity_check_group[0]) == 2 + + +def test_stim_text_to_pattern_maps_m_and_mpp_records_to_one_detector() -> None: + result = stim_text_to_pattern( + """ + MPP X0 + X_ERROR(0.01) 7 + M 7 + DETECTOR rec[-2] rec[-1] + """ + ) + + assert result.mpp_extractions[0].supports == (((0, "X"),),) + assert any( + isinstance(command, M) and isinstance(command.meas_basis, AxisMeasBasis) and command.meas_basis.axis == Axis.Z + for command in result.pattern.commands + ) + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + assert len(result.pattern.pauli_frame.parity_check_group[0]) == 2 + + +def test_stim_text_to_pattern_omits_noise_and_measurement_error_probabilities() -> None: + result = stim_text_to_pattern( + """ + DEPOLARIZE1(0.25) 0 + X_ERROR(0.125) 0 + MX(0.5) 0 + DETECTOR rec[-1] + """ + ) + + assert result.mpp_extractions == () + measurements = [command for command in result.pattern.commands if isinstance(command, M)] + assert len(measurements) == 1 + assert isinstance(measurements[0].meas_basis, AxisMeasBasis) + assert measurements[0].meas_basis.axis == Axis.X + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + + +def test_stim_text_to_pattern_preserves_ideal_herald_records_as_zero() -> None: + result = stim_text_to_pattern( + """ + MPP X0 + HERALDED_ERASE(0.25) 0 + DETECTOR rec[-2] rec[-1] + """ + ) + + assert result.mpp_extractions[0].supports == (((0, "X"),),) + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + assert len(result.pattern.pauli_frame.parity_check_group[0]) == 1 + + +def test_stim_text_to_pattern_preserves_cross_block_detector_records() -> None: + result = stim_text_to_pattern( + """ + MPP X0 + TICK + MPP Z0 + DETECTOR rec[-2] + """ + ) + + assert len(result.mpp_extractions) == 2 + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + assert len(result.pattern.pauli_frame.parity_check_group[0]) == 1 + + +def test_stim_text_to_pattern_tracks_all_record_types_with_global_indices() -> None: + result = stim_text_to_pattern( + """ + MPP X0 + HERALDED_ERASE(0.25) 2 + M 1 + TICK + MPP Z0 + DETECTOR rec[-4] rec[-3] rec[-2] rec[-1] + OBSERVABLE_INCLUDE(5) rec[-4] rec[-1] + """ + ) + + assert result.stim_to_qubit == {0: 0, 1: 1} + assert len(result.mpp_extractions) == 2 + for extraction in result.mpp_extractions: + assert extraction.detector_record_indices == (frozenset({0, 1, 2, 3}),) + assert extraction.logical_observable_record_indices == {5: frozenset({0, 3})} + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + assert len(result.pattern.pauli_frame.parity_check_group[0]) == 3 + assert len(result.pattern.pauli_frame.logical_observables[5]) == 2 + + +@pytest.mark.parametrize( + "text", + [ + "MPP X0\nDETECTOR rec[-1]", + "MPP X0\nTICK\nMPP X0\nDETECTOR rec[-1] rec[-2]", + ], +) +def test_stim_text_to_pattern_preserves_deterministic_mpp_detectors(text: str) -> None: + pattern = stim_text_to_pattern(text).pattern + compiled = stim.Circuit(stim_compile(pattern, emit_qubit_coords=False)) + + compiled.detector_error_model() + + +def test_stim_text_to_pattern_composes_mpp_output_into_next_mpp_input() -> None: + result = stim_text_to_pattern("MPP X0\nTICK\nMPP X0") + graph = result.pattern.pauli_frame.graphstate + + assert graph.number_of_nodes() == 7 + assert len(graph.meas_bases) == 6 + assert len(graph.input_node_indices) == 1 + assert len(graph.output_node_indices) == 1 + + +def test_stim_text_to_pattern_accepts_annotation_only_tick_block() -> None: + result = stim_text_to_pattern( + """ + MPP X0 + TICK + DETECTOR rec[-1] + """ + ) + + assert len(result.pattern.pauli_frame.parity_check_group) == 1 + + +@pytest.mark.parametrize("measurement", ["MPP X0", "M 0", "MX 0", "MY 0", "MXX 0 1"]) +def test_stim_text_to_pattern_rejects_mixed_measurement_and_unitary_block(measurement: str) -> None: + with pytest.raises(ValueError, match="separated from unitary gate instructions by TICK"): + stim_text_to_pattern(f"H 0\n{measurement}\n") + + +@pytest.mark.parametrize("instruction", ["R 0", "RX 0", "RY 0", "MR 0", "MRX 0", "MRY 0"]) +def test_stim_text_to_pattern_defers_reset_instructions(instruction: str) -> None: + with pytest.raises(ValueError, match="Unsupported Stim instruction"): + stim_text_to_pattern(instruction) + + +def test_stim_text_to_pattern_rejects_qubit_reuse_after_single_measurement() -> None: + with pytest.raises(ValueError, match="terminate those qubit lifetimes"): + stim_text_to_pattern("M 0\nTICK\nMPP X0") + + +def test_stim_text_to_pattern_rejects_unitary_reuse_after_single_measurement() -> None: + with pytest.raises(ValueError, match="terminate those qubit lifetimes"): + stim_text_to_pattern("M 0\nTICK\nH 0") + + +def test_stim_text_to_pattern_allows_disjoint_qubit_after_single_measurement() -> None: + result = stim_text_to_pattern("M 0\nTICK\nH 1") + measured_nodes = {command.node for command in result.pattern.commands if isinstance(command, M)} + simulator = PatternSimulator(result.pattern, SimulatorBackend.StateVector) + + simulator.simulate(rng=np.random.default_rng(3)) + + assert result.stim_to_qubit == {0: 0, 1: 1} + assert set(result.pattern.output_node_indices.values()) == {0, 1} + assert any(result.pattern.output_node_indices[node] == 0 for node in measured_nodes) + assert set(simulator.output_results) == {0} + assert np.isclose(abs(np.vdot([1, 0], simulator.state.state())), 1.0, atol=1e-9) + + +@pytest.mark.parametrize( + ("instruction", "expected_axis", "compiled_instruction"), + [ + ("MX !0", Axis.X, "MX !0"), + ("MY !0", Axis.Y, "MY !0"), + ("M !0", Axis.Z, "MZ !0"), + ], +) +def test_stim_text_to_pattern_preserves_inverted_single_measurement( + instruction: str, + expected_axis: Axis, + compiled_instruction: str, +) -> None: + result = stim_text_to_pattern(instruction) + measurements = [command for command in result.pattern.commands if isinstance(command, M)] + + assert result.mpp_extractions == () + assert len(measurements) == 1 + assert isinstance(measurements[0].meas_basis, AxisMeasBasis) + assert measurements[0].meas_basis.axis == expected_axis + assert measurements[0].meas_basis.sign == Sign.MINUS + assert compiled_instruction in stim_compile(result.pattern, emit_qubit_coords=False).splitlines() + + +def test_stim_text_to_pattern_rejects_inverted_pair_measurement_result() -> None: + with pytest.raises(ValueError, match="Signed MPP products are not supported"): + stim_text_to_pattern("MYY !0 1") + + +def test_stim_text_to_pattern_rejects_true_mpad_record() -> None: + with pytest.raises(ValueError, match="MPAD 1 records are not supported"): + stim_text_to_pattern("MPP X0\nMPAD 1") + + +def test_stim_text_to_pattern_rejects_record_before_beginning_of_time() -> None: + with pytest.raises(ValueError, match="before the beginning of time"): + stim_text_to_pattern("DETECTOR rec[-1]") diff --git a/tests/test_stim_mpp.py b/tests/test_stim_mpp.py index 68c65ab4..cf1d9316 100644 --- a/tests/test_stim_mpp.py +++ b/tests/test_stim_mpp.py @@ -50,7 +50,8 @@ def test_stabilizer_code_from_stim_mpp_preserves_3d_coordinates_when_requested() def test_stabilizer_code_from_stim_mpp_preserves_higher_dimensional_coordinates() -> None: extraction = stabilizer_code_from_stim_text( """ - QUBIT_COORDS(10, 20, 30, 40) 0 + QUBIT_COORDS(10, 20) 0 + QUBIT_COORDS(30, 40) 0 MPP X0 """, coord_dims=4, @@ -191,3 +192,8 @@ def test_stabilizer_code_from_stim_mpp_can_select_all_mpp_layers() -> None: assert extraction.supports == (((0, "X"),), ((0, "Z"),)) assert extraction.detector_rows == (frozenset({0, 1}),) assert extraction.detector_record_indices == (frozenset({0, 1}),) + + +def test_stabilizer_code_from_stim_mpp_rejects_record_before_beginning_of_time() -> None: + with pytest.raises(ValueError, match="before the beginning of time"): + stabilizer_code_from_stim_text("DETECTOR rec[-1]\nMPP X0")