Fix #177: Correct pattern export to OpenQASM 3

.github/workflows/typecheck.yml

@@ -12,7 +12,9 @@ concurrency:
   cancel-in-progress: true
 
 env:
-  python-version: "3.14"
+  # We stick to Python 3.13 since qiskit-aer is not available yet for Python 3.14.
+  # See https://github.com/Qiskit/qiskit-aer/issues/2378.
+  python-version: "3.13"
 
 jobs:
   mypy-pyright:

CHANGELOG.md

@@ -35,6 +35,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #394: The method `Circuit.transpile_measurements_to_z_axis` returns an equivalent circuit where all measurements are on the Z axis. This can be used to prepare a circuit for export to OpenQASM with `circuit_to_qasm3`.
 
+- #402: Support for Python 3.14.
+
 ### Fixed
 
 - #392: `Pattern.remove_input_nodes` is required before the `Pattern.perform_pauli_measurements` method to ensure input nodes are removed and fixed in the |+> state.
@@ -49,6 +51,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #389, #391: `Pattern.extract_opengraph` raises an exception if pattern has `N` commands which do not represent a |+> state.
 
+- #404: Fixed pattern export to OpenQASM 3. Compatibility with Qiskit
+  is ensured with normalization passed `incorporate_pauli_results` and
+  `single_qubit_domains`.
+
 - #409: Axis labels are shown when visualizing a pattern. Legend is placed outside the plot so that the graph remains visible.
 
 ### Changed
@@ -59,7 +65,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - #352, #394: Circuit measurements are now limited to axes X, Y, and Z.
 
-- #233, #398: The angle convention is now consistent across the
+- #233, #399: The angle convention is now consistent across the
   library: angles are represented as floats and expressed in units of
   π. In particular, angles that appear in parameters of circuit
   instructions are now expressed in units of π.

docs/source/optimization.rst

@@ -0,0 +1,19 @@
+Optimization passes
+===================
+
+:mod:`graphix.optimization` module
+++++++++++++++++++++++++++++++++++
+
+This module defines some optimization passes for patterns.
+
+.. currentmodule:: graphix.optimization
+
+.. autofunction:: standardize
+
+.. autoclass:: StandardizedPattern
+
+.. autofunction:: incorporate_pauli_results
+
+.. autofunction:: remove_useless_domains
+
+.. autofunction:: single_qubit_domains

docs/source/references.rst

@@ -16,3 +16,4 @@ Module reference
     channels
     random_objects
     open_graph
+    optimization

graphix/optimization.py

@@ -30,7 +30,7 @@
 from graphix.states import BasicStates
 
 if TYPE_CHECKING:
-    from collections.abc import Iterable, Mapping
+    from collections.abc import Callable, Iterable, Mapping
     from collections.abc import Set as AbstractSet
     from typing import Self
 
@@ -41,14 +41,15 @@ def standardize(pattern: Pattern) -> Pattern:
     """Return a standardized form to the given pattern.
 
     A standardized form is an equivalent pattern where the commands
-    appear in the following order: `N`, `E`, `M`, `Z`, `X`, `C`.
+    appear in the following order: ``N``, ``E``, ``M``, ``Z``, ``X``,
+    ``C``.
 
     Note that a standardized form does not always exist in presence of
-    `C` commands. For instance, there is no standardized form for the
+    ``C`` commands. For instance, there is no standardized form for the
     following pattern (written in the right-to-left convention):
-    `E(0, 1) C(0, H) N(1) N(0)`.
+    ``E(0, 1) C(0, H) N(1) N(0)``.
 
-    The function raises `NotImplementedError` if there is no
+    The function raises ``NotImplementedError`` if there is no
     standardized form. This behavior can change in the future.
 
 
@@ -61,6 +62,7 @@ def standardize(pattern: Pattern) -> Pattern:
     -------
     standardized : Pattern
         The standardized pattern, if it exists.
+
     """
     return StandardizedPattern.from_pattern(pattern).to_pattern()
 
@@ -97,7 +99,7 @@ class StandardizedPattern(_StandardizedPattern):
 
     Instances can be generated with the constructor from any
     compatible data structures, and an instance can be generated
-    directly from a pattern with the class method `from_pattern`.
+    directly from a pattern with the class method :meth:`from_pattern`.
 
     The constructor instantiates the ``Mapping`` fields as
     ``MappingProxyType`` objects over fresh dictionaries, ensuring
@@ -685,3 +687,37 @@ def remove_useless_domains(pattern: Pattern) -> Pattern:
             new_pattern.add(cmd)
     new_pattern.reorder_output_nodes(pattern.output_nodes)
     return new_pattern
+
+
+def single_qubit_domains(pattern: Pattern) -> Pattern:
+    """Return an equivalent pattern where domains contains at most one qubit."""
+    new_pattern = graphix.pattern.Pattern(input_nodes=pattern.input_nodes)
+    new_pattern.results = pattern.results
+
+    def decompose_domain(cmd: Callable[[int, set[int]], command.Command], node: int, domain: AbstractSet[int]) -> bool:
+        if len(domain) <= 1:
+            return False
+        for src in domain:
+            new_pattern.add(cmd(node, {src}))
+        return True
+
+    for cmd in pattern:
+        if cmd.kind == CommandKind.M:
+            replaced_s_domain = decompose_domain(command.X, cmd.node, cmd.s_domain)
+            replaced_t_domain = decompose_domain(command.Z, cmd.node, cmd.t_domain)
+            if replaced_s_domain or replaced_t_domain:
+                new_s_domain = set() if replaced_s_domain else cmd.s_domain
+                new_t_domain = set() if replaced_t_domain else cmd.t_domain
+                new_cmd = dataclasses.replace(cmd, s_domain=new_s_domain, t_domain=new_t_domain)
+                new_pattern.add(new_cmd)
+                continue
+        elif cmd.kind == CommandKind.X:
+            if decompose_domain(command.X, cmd.node, cmd.domain):
+                continue
+        elif cmd.kind == CommandKind.Z:
+            if decompose_domain(command.Z, cmd.node, cmd.domain):
+                continue
+        new_pattern.add(cmd)
+
+    new_pattern.reorder_output_nodes(pattern.output_nodes)
+    return new_pattern

graphix/pattern.py

@@ -27,6 +27,7 @@
 from graphix.measurements import Measurement, Outcome, PauliMeasurement, toggle_outcome
 from graphix.opengraph import OpenGraph
 from graphix.pretty_print import OutputFormat, pattern_to_str
+from graphix.qasm3_exporter import pattern_to_qasm3_lines
 from graphix.simulator import PatternSimulator
 from graphix.states import BasicStates
 from graphix.visualization import GraphVisualizer
@@ -41,6 +42,7 @@
     from graphix.flow.core import CausalFlow, GFlow
     from graphix.parameter import ExpressionOrSupportsFloat, Parameter
     from graphix.sim import Backend, BackendState, Data
+    from graphix.states import State
 
 
 _StateT_co = TypeVar("_StateT_co", bound="BackendState", covariant=True)
@@ -1498,27 +1500,21 @@ def draw_graph(
                 filename=filename,
             )
 
-    def to_qasm3(self, filename: Path | str) -> None:
+    def to_qasm3(self, filename: Path | str, input_state: dict[int, State] | State = BasicStates.PLUS) -> None:
         """Export measurement pattern to OpenQASM 3.0 file.
 
+        See :func:`graphix.qasm3_exporter.pattern_to_qasm3`.
+
         Parameters
         ----------
         filename : Path | str
-            file name to export to. example: "filename.qasm"
+            File name to export to. Example: ``"filename.qasm"``.
+
+        input_state : dict[int, State] | State, default BasicStates.PLUS
+            The initial state for each input node. Only ``|0⟩`` or ``|+⟩`` states are supported.
         """
         with Path(filename).with_suffix(".qasm").open("w", encoding="utf-8") as file:
-            file.write("// generated by graphix\n")
-            file.write("OPENQASM 3;\n")
-            file.write('include "stdgates.inc";\n')
-            file.write("\n")
-            if self.results != {}:
-                for i in self.results:
-                    res = self.results[i]
-                    file.write("// measurement result of qubit q" + str(i) + "\n")
-                    file.write("bit c" + str(i) + " = " + str(res) + ";\n")
-                    file.write("\n")
-            for cmd in self.__seq:
-                file.writelines(cmd_to_qasm3(cmd))
+            file.writelines(pattern_to_qasm3_lines(self, input_state=input_state))
 
     def is_parameterized(self) -> bool:
         """
@@ -1814,86 +1810,6 @@ def pauli_nodes(pattern: optimization.StandardizedPattern) -> tuple[list[tuple[c
     return pauli_node, non_pauli_node
 
 
-def cmd_to_qasm3(cmd: Command) -> Iterator[str]:
-    """Convert a command in the pattern into OpenQASM 3.0 statement.
-
-    Parameter
-    ---------
-    cmd : list
-        command [type:str, node:int, attr]
-
-    Yields
-    ------
-    string
-        translated pattern commands in OpenQASM 3.0 language
-
-    """
-    if cmd.kind == CommandKind.N:
-        qubit = cmd.node
-        yield "// prepare qubit q" + str(qubit) + "\n"
-        yield "qubit q" + str(qubit) + ";\n"
-        yield "h q" + str(qubit) + ";\n"
-        yield "\n"
-
-    elif cmd.kind == CommandKind.E:
-        qubits = cmd.nodes
-        yield "// entangle qubit q" + str(qubits[0]) + " and q" + str(qubits[1]) + "\n"
-        yield "cz q" + str(qubits[0]) + ", q" + str(qubits[1]) + ";\n"
-        yield "\n"
-
-    elif cmd.kind == CommandKind.M:
-        qubit = cmd.node
-        plane = cmd.plane
-        alpha = cmd.angle
-        sdomain = cmd.s_domain
-        tdomain = cmd.t_domain
-        yield "// measure qubit q" + str(qubit) + "\n"
-        yield "bit c" + str(qubit) + ";\n"
-        yield "float theta" + str(qubit) + " = 0;\n"
-        if plane == Plane.XY:
-            if sdomain:
-                yield "int s" + str(qubit) + " = 0;\n"
-                for sid in sdomain:
-                    yield "s" + str(qubit) + " += c" + str(sid) + ";\n"
-                yield "theta" + str(qubit) + " += (-1)**(s" + str(qubit) + " % 2) * (" + str(alpha) + " * pi);\n"
-            if tdomain:
-                yield "int t" + str(qubit) + " = 0;\n"
-                for tid in tdomain:
-                    yield "t" + str(qubit) + " += c" + str(tid) + ";\n"
-                yield "theta" + str(qubit) + " += t" + str(qubit) + " * pi;\n"
-            yield "p(-theta" + str(qubit) + ") q" + str(qubit) + ";\n"
-            yield "h q" + str(qubit) + ";\n"
-            yield "c" + str(qubit) + " = measure q" + str(qubit) + ";\n"
-            yield "h q" + str(qubit) + ";\n"
-            yield "p(theta" + str(qubit) + ") q" + str(qubit) + ";\n"
-            yield "\n"
-
-    # Use of == for mypy
-    elif cmd.kind == CommandKind.X or cmd.kind == CommandKind.Z:  # noqa: PLR1714
-        qubit = cmd.node
-        sdomain = cmd.domain
-        yield "// byproduct correction on qubit q" + str(qubit) + "\n"
-        yield "int s" + str(qubit) + " = 0;\n"
-        for sid in sdomain:
-            yield "s" + str(qubit) + " += c" + str(sid) + ";\n"
-        yield "if(s" + str(qubit) + " % 2 == 1){\n"
-        if cmd.kind == CommandKind.X:
-            yield "\t x q" + str(qubit) + ";\n}\n"
-        else:
-            yield "\t z q" + str(qubit) + ";\n}\n"
-        yield "\n"
-
-    elif cmd.kind == CommandKind.C:
-        qubit = cmd.node
-        yield "// Clifford operations on qubit q" + str(qubit) + "\n"
-        for op in cmd.clifford.qasm3:
-            yield str(op) + " q" + str(qubit) + ";\n"
-        yield "\n"
-
-    else:
-        raise ValueError(f"invalid command {cmd}")
-
-
 def assert_permutation(original: list[int], user: list[int]) -> None:
     """Check that the provided `user` node list is a permutation from `original`."""
     node_set = set(user)