Random generator of open graphs with gflow

graphix/random_objects.py

@@ -5,18 +5,24 @@
 import functools
 from typing import TYPE_CHECKING
 
+import networkx as nx
 import numpy as np
 import numpy.typing as npt
 import scipy.linalg
 from scipy.stats import unitary_group
 
+from graphix._linalg import MatGF2
 from graphix.channels import KrausChannel, KrausData
+from graphix.fundamentals import Plane
+from graphix.measurements import Measurement
+from graphix.opengraph import OpenGraph
 from graphix.ops import Ops
+from graphix.parameter import Placeholder
 from graphix.rng import ensure_rng
 from graphix.transpiler import Circuit
 
 if TYPE_CHECKING:
-    from collections.abc import Iterable, Iterator
+    from collections.abc import Iterable, Iterator, Mapping
 
     from numpy.random import Generator
 
@@ -423,3 +429,139 @@ def rand_circuit(
             ind = rng.integers(len(gate_choice))
             gate_choice[ind](j)
     return circuit
+
+
+def rand_og_gflow(n: int, n_o: int, rng: Generator, meas_planes: Mapping[int, Plane] | None = None) -> OpenGraph:
+    """Return random open graph with gflow.
+
+    Parameters
+    ----------
+    n : int
+        Number of nodes in open graph.
+    n_o : int
+        Number of output nodes in open graph
+    rng : numpy.random.Generator
+        Random number generator.
+    meas_planes : Mapping[int, Plane] | None
+        Measurement planes (`value`) associated with each non-output node (`key`). If `None`, they are chosen at random. Optional (defaults to `None`).
+
+    Returns
+    -------
+    graphix.opengraph.OpenGraph
+        Random open graph with gflow.
+
+    Notes
+    -----
+    This function implements the PGA algorithm presented in [1].
+
+    References
+    ----------
+    [1] "Generation of random open graph with gflow" Rajarsi Pal, Harold Ollivier, Graphix Workshop, Inria, 2024.
+    """
+    n_no = n - n_o  # Number of non-output nodes.
+    if meas_planes is not None and len(meas_planes) != n_no:
+        raise ValueError(f"Number of items in `meas_planes` should be equal to number of non-outputs, n - no = {n_no}")
+    if meas_planes is None:
+        plane_choice = [Plane.XY, Plane.XZ, Plane.YZ]
+        meas_planes = {i: plane_choice[rng.integers(3)] for i in range(n_no)}
+
+    shift = max(meas_planes) + 1
+    idx_node_mapping = {i: i + shift for i in range(n_o)}  # Index to node mapping in the adjacency matrix.
+
+    c_matrix = np.zeros((n_no + 1, n), dtype=np.uint8).view(MatGF2)
+    adj_matrix = np.zeros((n, n), dtype=np.uint8).view(MatGF2)
+
+    for k, (node, plane) in enumerate(meas_planes.items()):
+        idx_node_mapping[k + n_o] = node
+        c_matrix, adj_matrix = _add_vertex(c_matrix, adj_matrix, k, plane, n_o, rng)
+    graph = nx.from_numpy_array(adj_matrix)
+    graph = nx.relabel_nodes(graph, idx_node_mapping)
+    measurements = {i: Measurement(Placeholder("Angle"), plane) for i, plane in meas_planes.items()}
+    outputs = [idx_node_mapping[i] for i in range(n_o)]
+    inputs = [idx_node_mapping[i] for i, col in enumerate(c_matrix.T) if not np.any(col)]
+
+    return OpenGraph(inside=graph, measurements=measurements, inputs=inputs, outputs=outputs)
+
+
+def _add_vertex(
+    c_matrix: MatGF2, adj_matrix: MatGF2, k: int, plane: Plane, n_o: int, rng: Generator
+) -> tuple[MatGF2, MatGF2]:
+    """Add vertex to open graph.
+
+    Parameters
+    ----------
+    c_matrix : graphix._linalg.MatGF2
+        Correction matrix.
+    adj_matrix: graphix._linalg.MatGF2
+        Graph adjacency matrix.
+    k : int
+        Index of the updated row (and column) of the adjacency matrix.
+    plane : graphix.fundamentals.Plane
+        Measurement plane of the node under evaluation.
+    n_o: int
+        Number of output nodes in open graph.
+    rng : numpy.random.Generator
+        Random number generator.
+
+    Returns
+    -------
+    c_matrix : MatGF2
+        Updated correction matrix.
+    adj_matrix: MatGF2
+        Updated graph adjacency matrix.
+    """
+    k_shift = k + n_o
+    kernel = c_matrix[: k + 1, :k_shift].view(MatGF2).null_space()
+    # We pick a random linear combination of kernel-basis vector (with random coefficients)
+    mask = rng.choice([True, False], size=kernel.shape[0])
+    # Ensure that we have at least one non-zero coefficient.
+    if not np.any(mask):
+        mask[rng.integers(len(mask))] = True
+    g_vect = np.bitwise_xor.reduce(kernel[mask], axis=0)  # `g_vect` has shape (k_shift, )
+    c_vect = np.empty(k_shift + 1, dtype=np.uint8)  # `c_vect` has shape (k_shift + 1 , )
+
+    if plane == Plane.XY:
+        c_dot, c_k = 1, 0
+    elif plane == Plane.YZ:
+        c_dot, c_k = 0, 1
+    elif plane == Plane.XZ:
+        c_dot, c_k = 1, 1
+
+    c_vect[:k_shift] = _generate_rnd_gf2_vec(g_vect, c_dot, rng)
+    c_vect[-1] = c_k
+
+    adj_matrix[k_shift, :k_shift] = g_vect
+    adj_matrix[:k_shift, k_shift] = g_vect
+    c_matrix[k + 1, : k_shift + 1] = c_vect
+
+    return c_matrix.view(MatGF2), adj_matrix.view(MatGF2)
+
+
+def _generate_rnd_gf2_vec(g: npt.NDArray[np.uint8], c: int, rng: Generator) -> npt.NDArray[np.uint8]:
+    r"""Generate a random vector :math:`x` over :math:`\mathbb F_2` such that :math:`x \dot g = c`.
+
+    Parameters
+    ----------
+    g : npt.NDArray[np.uint8]
+    c : int (0 or 1)
+    rng : Generator
+
+    Returns
+    -------
+    npt.NDArray[np.uint8]
+        A random vector with `len(g)` elements and the appropriate constraints.
+    """
+    result = rng.integers(0, 2, size=len(g), dtype=np.uint8)
+
+    idx_nonzero = np.flatnonzero(g)
+
+    if idx_nonzero.size == 0:
+        if c != 0:
+            raise ValueError(r"It does not exist an `x` such that  $x \cdot g = c$.")
+        return result
+
+    # If the random vector `result` does not fulfil the constraint, it suffices to flip bit `k`, such that `g[k] != 0`.
+    if np.bitwise_xor.reduce(result[idx_nonzero]) != c:
+        result[rng.choice(idx_nonzero)] ^= 1
+
+    return result

tests/test_random_utilities.py

@@ -8,6 +8,8 @@
 import graphix.random_objects as randobj
 from graphix import linalg_validations as lv
 from graphix.channels import KrausChannel
+from graphix.find_pflow import find_pflow
+from graphix.fundamentals import Plane
 from graphix.ops import Ops
 from graphix.sim.density_matrix import DensityMatrix
 
@@ -179,3 +181,49 @@ def test_random_pauli_channel_fail(self, fx_rng: Generator) -> None:
 
         with pytest.raises(ValueError):
             randobj.rand_pauli_channel_kraus(dim=2**nqb + 1, rank=rk, rng=fx_rng)
+
+
+class TestRandomOpenGraph:
+    """A class to group all the tests related to graphix.random_objects.rand_og_flow."""
+
+    def test_generate_rnd_gf2_vec(self, fx_rng: Generator) -> None:
+        """Test that returned vector fulfils constraints and are uniformly disrtibuted."""
+        n_shots = 1000  # number of shots
+        g = np.array([1, 0, 0], dtype=np.uint8)
+        c = 1
+
+        x_ref = {
+            0: np.array([1, 1, 0], dtype=np.uint8),
+            1: np.array([1, 1, 1], dtype=np.uint8),
+            2: np.array([1, 0, 0], dtype=np.uint8),
+            3: np.array([1, 0, 1], dtype=np.uint8),
+        }
+        counts = {0: 0, 1: 0, 2: 0, 3: 0}
+
+        for _ in range(n_shots):
+            x = randobj._generate_rnd_gf2_vec(g, c, fx_rng)
+            for i in range(4):
+                if np.all(x == x_ref[i]):
+                    counts[i] += 1
+                    break
+            else:
+                pytest.fail(f"Generated impossible random vector x = {x}")
+
+        for count in counts.values():
+            assert count / n_shots == pytest.approx(0.25, abs=1 / np.sqrt(n_shots))
+
+    def test_rand_og_gflow(self, fx_rng: Generator) -> None:
+        meas_planes = {0: Plane.XY, 1: Plane.XZ, 3: Plane.YZ}
+        n = 5
+        n_o = 2
+        og = randobj.rand_og_gflow(n, n_o, fx_rng, meas_planes)
+        assert og.inside.number_of_nodes() == n
+        assert og.outputs == [4, 5]
+        assert {node: meas.plane for node, meas in og.measurements.items()} == meas_planes
+        assert find_pflow(og) is not None
+
+        n = 10
+        n_o = 4
+        for _ in range(10):
+            og = randobj.rand_og_gflow(n, n_o, fx_rng)
+            assert find_pflow(og) is not None