Add ZX-calculus backend

docs/source/gflow_utils.rst

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

docs/source/graphstate.rst

@@ -24,11 +24,3 @@ Functions
 .. autofunction:: graphqomb.graphstate.compose
 .. autofunction:: graphqomb.graphstate.bipartite_edges
 .. autofunction:: graphqomb.graphstate.odd_neighbors
-
-Auxiliary Classes
-------------------
-.. autoclass:: graphqomb.graphstate.LocalCliffordExpansion
-    :members:
-
-.. autoclass:: graphqomb.graphstate.ExpansionMaps
-    :members:

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,32 @@
+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.complete_graph_edges
+
+Auxiliary Classes
+------------------
+
+.. autoclass:: graphqomb.zxgraphstate.InputLocalCliffordExpansion
+   :members:
+
+.. autoclass:: graphqomb.zxgraphstate.OutputLocalCliffordExpansion
+   :members:
+
+.. autoclass:: graphqomb.zxgraphstate.ExpansionMaps
+   :members:

examples/zxgraph_simplification.py

@@ -0,0 +1,136 @@
+"""
+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
+
+# %%
+# Prepare an initial random graph state with flow
+graph, flow = generate_random_flow_graph(width=3, depth=4, edge_p=0.5)
+zx_graph, _ = ZXGraphState.from_base_graph_state(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)
+
+
+# %%
+# 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()
+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.
+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/euler.py

@@ -213,41 +213,70 @@ def conjugate(self) -> LocalClifford:
         return LocalClifford(-self.gamma, -self.beta, -self.alpha)
 
 
-def meas_basis_info(vector: NDArray[np.complex128]) -> tuple[Plane, float]:
+def _meas_basis_candidates(vector: NDArray[np.complex128]) -> list[tuple[Plane, float]]:
+    r"""Return candidate measurement planes and angles corresponding to a vector.
+
+    Parameters
+    ----------
+    vector : `numpy.typing.NDArray`\[`numpy.complex128`\]
+        1 qubit state vector
+
+    Returns
+    -------
+    `list`\[`tuple`\[`Plane`, `float`\]\]
+        candidate measurement planes and angles
+    """
+    theta, phi = bloch_sphere_coordinates(vector)
+    candidates: list[tuple[Plane, float]] = []
+    if is_clifford_angle(phi):
+        # YZ or XZ plane
+        if is_close_angle(2 * phi, 0):  # 0 or pi
+            xz_theta = -theta if is_close_angle(phi, math.pi) else theta
+            candidates.append((Plane.XZ, xz_theta))
+            if is_close_angle(phi, 0) and is_close_angle(2 * theta, 0):
+                candidates.append((Plane.YZ, theta))
+        if is_close_angle(2 * (phi - math.pi / 2), 0):
+            yz_theta = -theta if is_close_angle(phi, 3 * math.pi / 2) else theta
+            candidates.append((Plane.YZ, yz_theta))
+    if is_clifford_angle(theta) and not is_clifford_angle(theta / 2):
+        # XY plane
+        phi = phi + math.pi if is_close_angle(theta, 3 * math.pi / 2) else phi
+        candidates.append((Plane.XY, phi))
+
+    return candidates
+
+
+def meas_basis_info(vector: NDArray[np.complex128], expected_plane: Plane | None = None) -> tuple[Plane, float]:
     r"""Return the measurement plane and angle corresponding to a vector.
 
     Parameters
     ----------
     vector : `numpy.typing.NDArray`\[`numpy.complex128`\]
         1 qubit state vector
+    expected_plane : `Plane` | `None`, optional
+        expected measurement plane to preserve gflow existence, by default None
 
     Returns
     -------
-    `tuple`\[`Plane`, `float`]
+    `tuple`\[`Plane`, `float`\]
         measurement plane and angle
 
     Raises
     ------
     ValueError
         if the vector does not lie on any of 3 planes
     """
-    theta, phi = bloch_sphere_coordinates(vector)
-    if is_clifford_angle(phi):
-        # YZ or XZ plane
-        if is_clifford_angle(phi / 2):  # 0 or pi
-            if is_close_angle(phi, math.pi):
-                theta = -theta
-            return Plane.XZ, theta
-        if is_close_angle(phi, 3 * math.pi / 2):
-            theta = -theta
-        return Plane.YZ, theta
-    if is_clifford_angle(theta) and not is_clifford_angle(theta / 2):
-        # XY plane
-        if is_close_angle(theta, 3 * math.pi / 2):
-            phi += math.pi
-        return Plane.XY, phi
-    msg = "The vector does not lie on any of 3 planes"
-    raise ValueError(msg)
+    candidates = _meas_basis_candidates(vector)
+
+    if not candidates:
+        msg = "The vector does not lie on any of 3 planes"
+        raise ValueError(msg)
+
+    if expected_plane is not None:
+        for plane, angle in candidates:
+            if plane == expected_plane:
+                return plane, angle
+    return candidates[0]
 
 
 # TODO(masa10-f): Algebraic backend for this computation(#023)
@@ -275,7 +304,7 @@ def update_lc_lc(lc1: LocalClifford, lc2: LocalClifford) -> LocalClifford:
 
 
 # TODO(masa10-f): Algebraic backend for this computation(#023)
-def update_lc_basis(lc: LocalClifford, basis: MeasBasis) -> PlannerMeasBasis:
+def update_lc_basis(lc: LocalClifford, basis: MeasBasis, expected_plane: Plane | None = None) -> PlannerMeasBasis:
     """Update a `MeasBasis` object with an action of `LocalClifford` object.
 
     Parameters
@@ -284,17 +313,19 @@ def update_lc_basis(lc: LocalClifford, basis: MeasBasis) -> PlannerMeasBasis:
         `LocalClifford`
     basis : `MeasBasis`
         `MeasBasis`
+    expected_plane : `Plane` | `None`, optional
+        expected measurement plane to preserve gflow existence, by default None
 
     Returns
     -------
     `PlannerMeasBasis`
         updated `PlannerMeasBasis`
     """
-    matrix = lc.matrix()
+    matrix = lc.matrix().conjugate().T
     vector = basis.vector()
 
     updated_vector = np.asarray(matrix @ vector, dtype=np.complex128)
-    plane, angle = meas_basis_info(updated_vector)
+    plane, angle = meas_basis_info(updated_vector, expected_plane=expected_plane)
     return PlannerMeasBasis(plane, angle)
 
 

graphqomb/gflow_utils.py

@@ -0,0 +1,68 @@
+"""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.
+"""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import networkx as nx
+from swiflow import gflow
+from swiflow.common import Plane as SfPlane
+from typing_extensions import assert_never
+
+from graphqomb.common import Plane
+
+if TYPE_CHECKING:
+    from networkx import Graph as NxGraph
+
+    from graphqomb.graphstate import BaseGraphState
+
+
+def gflow_wrapper(graphstate: BaseGraphState) -> dict[int, set[int]]:
+    """Utilize ``swiflow.gflow`` to search gflow.
+
+    Parameters
+    ----------
+    graphstate : `BaseGraphState`
+        graph state to find gflow
+
+    Returns
+    -------
+    ``dict[int, set[int]]``
+        gflow object
+
+    Raises
+    ------
+    ValueError
+        If no gflow is found
+    """
+    graph: NxGraph[int] = 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:
+            assert_never(plane)
+
+    gflow_object = gflow.find(
+        graph, graphstate.input_node_indices.keys(), graphstate.output_node_indices.keys(), swiflow_planes
+    )
+    if gflow_object is None:
+        msg = "No flow found"
+        raise ValueError(msg)
+
+    return gflow_object.f