Add ZX-calculus backend

docs/source/gflow_utils.rst

@@ -0,0 +1,17 @@
+Gflow Utils
+===========
+
+:mod:`graphqomb.gflow_utils` module
++++++++++++++++++++++++++++++++++++
+
+.. automodule:: graphqomb.gflow_utils
+
+Functions
+---------
+
+.. autofunction:: graphqomb.gflow_utils.gflow_wrapper
+
+Constants
+---------
+
+.. autodata:: graphqomb.gflow_utils._EQUIV_MEAS_BASIS_MAP

docs/source/references.rst

@@ -15,10 +15,12 @@ Module reference
     random_objects
     feedforward
     focus_flow
+    gflow_utils
     command
     pattern
     pauli_frame
     qompiler
     scheduler
     stim_compiler
     visualizer
+    zxgraphstate

docs/source/zxgraphstate.rst

@@ -0,0 +1,21 @@
+ZXGraphState
+============
+
+:mod:`graphqomb.zxgraphstate` module
++++++++++++++++++++++++++++++++++++++
+
+.. automodule:: graphqomb.zxgraphstate
+
+ZX Graph State
+--------------
+
+.. autoclass:: graphqomb.zxgraphstate.ZXGraphState
+   :members:
+   :show-inheritance:
+   :member-order: bysource
+
+Functions
+---------
+
+.. autofunction:: graphqomb.zxgraphstate.to_zx_graphstate
+.. autofunction:: graphqomb.zxgraphstate.complete_graph_edges

examples/zxgraph_simplification.py

@@ -0,0 +1,154 @@
+"""
+Basic example of simplifying a ZX-diagram.
+==========================================
+
+By using the full_reduce method,
+we can remove all the internal Clifford nodes and some non-Clifford nodes from the graph state,
+which generates a simpler ZX-diagram.
+This example is a simple demonstration of the simplification process.
+
+Note that as a result of the simplification, local Clifford operations are applied to the input/output nodes.
+"""
+
+# %%
+
+from __future__ import annotations
+
+from copy import deepcopy
+
+import numpy as np
+
+from graphqomb.gflow_utils import gflow_wrapper
+from graphqomb.qompiler import qompile
+from graphqomb.random_objects import generate_random_flow_graph
+from graphqomb.simulator import PatternSimulator, SimulatorBackend
+from graphqomb.visualizer import visualize
+from graphqomb.zxgraphstate import ZXGraphState, to_zx_graphstate
+
+FlowLike = dict[int, set[int]]
+
+# %%
+# Prepare an initial random graph state with flow
+graph, flow = generate_random_flow_graph(width=3, depth=4, edge_p=0.5)
+zx_graph, _ = to_zx_graphstate(graph)
+visualize(zx_graph)
+
+# %%
+# We can compile the graph state into a measurement pattern, simulate it, and get the resulting statevector.
+pattern = qompile(zx_graph, flow)
+sim = PatternSimulator(pattern, backend=SimulatorBackend.StateVector)
+sim.simulate()
+statevec_original = sim.state
+
+
+# %%
+def print_boundary_lcs(zxgraph: ZXGraphState) -> None:
+    lc_map = zxgraph.local_cliffords
+    for node in zxgraph.input_node_indices | zxgraph.output_node_indices:
+        # check lc on input and output nodes
+        lc = lc_map.get(node, None)
+        if lc is not None:
+            if node in zxgraph.input_node_indices:
+                print(f"Input node {node} has local Clifford: alpha={lc.alpha}, beta={lc.beta}, gamma={lc.gamma}")
+            else:
+                print(f"Output node {node} has local Clifford: alpha={lc.alpha}, beta={lc.beta}, gamma={lc.gamma}")
+        else:
+            print(f"Node {node} has no local Clifford.")
+
+
+def print_meas_bses(graph: ZXGraphState) -> None:
+    print("node | plane | angle (/pi)")
+    for node in graph.input_node_indices:
+        print(f"{node} (input)", graph.meas_bases[node].plane, graph.meas_bases[node].angle / np.pi)
+    for node in graph.physical_nodes - set(graph.input_node_indices) - set(graph.output_node_indices):
+        print(node, graph.meas_bases[node].plane, graph.meas_bases[node].angle / np.pi)
+    for node in graph.output_node_indices:
+        print(f"{node} (output)", "-", "-")
+
+
+# %%
+print_boundary_lcs(zx_graph)
+
+# %%
+# Initial graph state before simplification
+print_meas_bses(zx_graph)
+
+
+# %%
+# Simplify the graph state by full_reduce method
+zx_graph_smp = deepcopy(zx_graph)
+zx_graph_smp.full_reduce()
+
+# %%
+# Simplified graph state after full_reduce.
+visualize(zx_graph_smp)
+print_meas_bses(zx_graph_smp)
+print_boundary_lcs(zx_graph_smp)
+
+# %%
+# NOTE:
+# At first glance, the input/output nodes appear to remain unaffected.
+# However, note that a local Clifford operation is actually applied as a result of the action of the full_reduce method.
+
+# If you visualize the graph state after executing the `expand_local_cliffords` method,
+# you will see additional nodes connected to the former input/output nodes.
+
+
+# %%
+# Let us compare the graph state before and after simplification.
+# We simulate the pattern obtained from the simplified graph state.
+# Note that we need to call the `expand_local_cliffords` method before generating the pattern to get the gflow.
+
+zx_graph_smp.expand_local_cliffords()
+zx_graph_smp.to_xy()  # to improve gflow search performance
+zx_graph_smp.to_xz()  # to improve gflow search performance
+print("input_node_indices: ", set(zx_graph_smp.input_node_indices))
+print("output_node_indices: ", set(zx_graph_smp.output_node_indices))
+print("local_cliffords: ", zx_graph_smp.local_cliffords)
+
+print_meas_bses(zx_graph_smp)
+visualize(zx_graph_smp)
+print_boundary_lcs(zx_graph_smp)
+
+# %%
+# Now we can obtain the gflow for the simplified graph state.
+# Then, we compile the simplified graph state into a measurement pattern,
+# simulate it, and get the resulting statevector.
+
+# NOTE:
+# gflow_wrapper does not support graph states with multiple subgraph structures in the gflow search wrapper below.
+# Hence, in case you fail, ensure that the simplified graph state consists of a single connected component.
+# To calculate the graph states with multiple subgraph structures,
+# you need to calculate gflow for each connected component separately.
+gflow_smp = gflow_wrapper(zx_graph_smp)
+pattern_smp = qompile(zx_graph_smp, gflow_smp)
+sim_smp = PatternSimulator(pattern_smp, backend=SimulatorBackend.StateVector)
+sim_smp.simulate()
+
+# %%
+statevec_smp = sim_smp.state
+# %%
+# normalization check
+print("norm of original statevector:", np.linalg.norm(statevec_original.state()))
+print("norm of simplified statevector:", np.linalg.norm(statevec_smp.state()))
+
+# %%
+# Finally, we compare the expectation values of random observables before and after simplification.
+rng = np.random.default_rng()
+for i in range(len(zx_graph.input_node_indices)):
+    rand_mat = rng.random((2, 2)) + 1j * rng.random((2, 2))
+    rand_mat += rand_mat.T.conj()
+    exp = statevec_original.expectation(rand_mat, [i])
+    exp_cr = statevec_smp.expectation(rand_mat, [i])
+    print("Expectation values for rand_mat\n===============================")
+    print("rand_mat: \n", rand_mat)
+    print("Original: \t\t", exp)
+    print("After simplification: \t", exp_cr)
+
+print("norm: ", np.linalg.norm(statevec_original.state()), np.linalg.norm(statevec_smp.state()))
+print("data shape: ", statevec_original.state().shape, statevec_smp.state().shape)
+psi_org = statevec_original.state()
+psi_smp = statevec_smp.state()
+print("inner product: ", np.abs(np.vdot(psi_org, psi_smp)))
+
+# %%

graphqomb/common.py

@@ -397,3 +397,45 @@ def meas_basis(plane: Plane, angle: float) -> NDArray[np.complex128]:
     else:
         typing_extensions.assert_never(plane)
     return basis.astype(np.complex128)
+
+
+def basis2tuple(meas_basis: MeasBasis) -> tuple[Plane, float]:
+    r"""Return the key (tuple[Plane, float]) for _EQUIV_MEAS_BASIS_MAP from a measurement basis.
+
+    Parameters
+    ----------
+    meas_basis : `MeasBasis`
+        measurement basis
+
+    Returns
+    -------
+    tuple[Plane, float]
+        key for _EQUIV_MEAS_BASIS_MAP
+    """
+    angle = round_clifford_angle(meas_basis.angle)
+    return (meas_basis.plane, angle)
+
+
+def round_clifford_angle(angle: float, atol: float = 1e-9) -> float:
+    r"""Round the Clifford angle numerically to the nearest Clifford angle.
+
+    Parameters
+    ----------
+    angle : `float`
+        angle in radians
+    atol : `float`, optional
+        absolute tolerance, by default 1e-9
+
+    Returns
+    -------
+    angle : `float`
+        For Clifford angles, the rounded angle.
+        If the angle is not a Clifford angle, return the original angle.
+    """
+    clifford_angles = [0.0, np.pi / 2, np.pi, 3 * np.pi / 2]
+    for ca in clifford_angles:
+        if is_close_angle(angle, ca, atol):
+            angle = ca
+            break
+
+    return angle

graphqomb/euler.py

@@ -290,7 +290,7 @@ def update_lc_basis(lc: LocalClifford, basis: MeasBasis) -> PlannerMeasBasis:
     `PlannerMeasBasis`
         updated `PlannerMeasBasis`
     """
-    matrix = lc.matrix()
+    matrix = lc.matrix().conjugate().T
     vector = basis.vector()
 
     updated_vector = np.asarray(matrix @ vector, dtype=np.complex128)

graphqomb/gflow_utils.py

@@ -0,0 +1,111 @@
+"""Utilities for generalized flow (gflow) computation.
+
+This module provides:
+
+- `gflow_wrapper`: Thin adapter around ``swiflow.gflow`` so that gflow can be computed directly
+    from a `BaseGraphState` instance.
+- `_EQUIV_MEAS_BASIS_MAP`: A mapping between equivalent measurement bases used to improve gflow finding performance.
+"""
+
+from __future__ import annotations
+
+import math
+from typing import TYPE_CHECKING
+
+import networkx as nx
+from swiflow import gflow
+from swiflow.common import Plane as SfPlane
+
+from graphqomb.common import Plane, PlannerMeasBasis
+
+if TYPE_CHECKING:
+    from typing import Any
+
+    from networkx import Graph as NxGraph
+
+    from graphqomb.graphstate import BaseGraphState
+
+    FlowLike = dict[int, set[int]]
+
+
+def gflow_wrapper(graphstate: BaseGraphState) -> FlowLike:
+    """Utilize ``swiflow.gflow`` to search gflow.
+
+    Parameters
+    ----------
+    graphstate : `BaseGraphState`
+        graph state to find gflow
+
+    Returns
+    -------
+    ``FlowLike``
+        gflow object
+
+    Raises
+    ------
+    ValueError
+        If no gflow is found
+
+    Notes
+    -----
+    This wrapper does not support graph states with multiple subgraph structures.
+    """
+    graph: NxGraph[Any] = nx.Graph()
+    graph.add_nodes_from(graphstate.physical_nodes)
+    graph.add_edges_from(graphstate.physical_edges)
+
+    bases = graphstate.meas_bases
+    planes = {node: bases[node].plane for node in bases}
+    swiflow_planes: dict[int, SfPlane] = {}
+    for node, plane in planes.items():
+        if plane == Plane.XY:
+            swiflow_planes[node] = SfPlane.XY
+        elif plane == Plane.YZ:
+            swiflow_planes[node] = SfPlane.YZ
+        elif plane == Plane.XZ:
+            swiflow_planes[node] = SfPlane.XZ
+        else:
+            msg = f"No match {plane}"
+            raise ValueError(msg)
+
+    gflow_object = gflow.find(
+        graph, set(graphstate.input_node_indices), set(graphstate.output_node_indices), swiflow_planes
+    )
+    if gflow_object is None:
+        msg = "No flow found"
+        raise ValueError(msg)
+
+    gflow_obj = gflow_object.f
+
+    return {node: {child for child in children if child != node} for node, children in gflow_obj.items()}
+
+
+#: Mapping between equivalent measurement bases.
+#:
+#: This map is used to replace a measurement basis by an equivalent one
+#: to improve gflow search performance.
+#:
+#: Key:
+#:   ``(Plane, angle)`` where angle is in radians.
+#: Value:
+#:   :class:`~graphqomb.common.PlannerMeasBasis`.
+_EQUIV_MEAS_BASIS_MAP: dict[tuple[Plane, float], PlannerMeasBasis] = {
+    # (XY, 0) <-> (XZ, pi/2)
+    (Plane.XY, 0.0): PlannerMeasBasis(Plane.XZ, 0.5 * math.pi),
+    (Plane.XZ, 0.5 * math.pi): PlannerMeasBasis(Plane.XY, 0.0),
+    # (XY, pi/2) <-> (YZ, pi/2)
+    (Plane.XY, 0.5 * math.pi): PlannerMeasBasis(Plane.YZ, 0.5 * math.pi),
+    (Plane.YZ, 0.5 * math.pi): PlannerMeasBasis(Plane.XY, 0.5 * math.pi),
+    # (XY, -pi/2) == (XY, 3pi/2) <-> (YZ, 3pi/2)
+    (Plane.XY, 1.5 * math.pi): PlannerMeasBasis(Plane.YZ, 1.5 * math.pi),
+    (Plane.YZ, 1.5 * math.pi): PlannerMeasBasis(Plane.XY, 1.5 * math.pi),
+    # (XY, pi) <-> (XZ, -pi/2) == (XZ, 3pi/2)
+    (Plane.XY, math.pi): PlannerMeasBasis(Plane.XZ, 1.5 * math.pi),
+    (Plane.XZ, 1.5 * math.pi): PlannerMeasBasis(Plane.XY, math.pi),
+    # (XZ, 0) <-> (YZ, 0)
+    (Plane.XZ, 0.0): PlannerMeasBasis(Plane.YZ, 0.0),
+    (Plane.YZ, 0.0): PlannerMeasBasis(Plane.XZ, 0.0),
+    # (XZ, pi) <-> (YZ, pi)
+    (Plane.XZ, math.pi): PlannerMeasBasis(Plane.YZ, math.pi),
+    (Plane.YZ, math.pi): PlannerMeasBasis(Plane.XZ, math.pi),
+}