diff --git a/graphix/random_objects.py b/graphix/random_objects.py index 5249c832..19d06f24 100644 --- a/graphix/random_objects.py +++ b/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 diff --git a/tests/test_random_utilities.py b/tests/test_random_utilities.py index 539874fe..94059807 100644 --- a/tests/test_random_utilities.py +++ b/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