diff --git a/CHANGELOG.md b/CHANGELOG.md index 2df08aadb..763cc30b0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- #393 + - Introduced new method `graphix.optimization.StandardizedPattern.extract_partial_order_layer` which constructs a partial order layering from the dependency domains of M, X and Z commands. + - Introduced new methods `graphix.optimization.StandardizedPattern.extract_causal_flow`, `graphix.optimization.StandardizedPattern.extract_gflow` which respectively attempt to extract a causal flow and a gflow from a standardized pattern. + - Introduced new wrapper methods in the `Pattern` class: `graphix.pattern.Pattern.extract_partial_order_layers`, `graphix.pattern.Pattern.extract_causal_flow` and `graphix.pattern.Pattern.extract_gflow`. + - Introduced new module `graphix.flow._partial_order` with the function :func:`compute_topological_generations`. + - #385 - Introduced `graphix.flow.core.XZCorrections.check_well_formed` which verifies the correctness of an XZ-corrections instance and raises an exception if incorrect. - Added XZ-correction exceptions to module `graphix.flow.core.exceptions`. diff --git a/graphix/flow/_find_gpflow.py b/graphix/flow/_find_gpflow.py index 65e3ac917..ea903a1cd 100644 --- a/graphix/flow/_find_gpflow.py +++ b/graphix/flow/_find_gpflow.py @@ -20,6 +20,7 @@ from typing_extensions import override from graphix._linalg import MatGF2, solve_f2_linear_system +from graphix.flow._partial_order import compute_topological_generations from graphix.fundamentals import AbstractMeasurement, AbstractPlanarMeasurement, Axis, Plane from graphix.sim.base_backend import NodeIndex @@ -533,7 +534,7 @@ def _compute_correction_matrix_general_case( return c_prime_matrix.mat_mul(cb_matrix) -def _compute_topological_generations(ordering_matrix: MatGF2) -> list[list[int]] | None: +def _try_ordering_matrix_to_topological_generations(ordering_matrix: MatGF2) -> tuple[frozenset[int], ...] | None: """Stratify the directed acyclic graph (DAG) represented by the ordering matrix into generations. Parameters @@ -543,8 +544,8 @@ def _compute_topological_generations(ordering_matrix: MatGF2) -> list[list[int]] Returns ------- - list[list[int]] - topological generations. Integers represent the indices of the matrix `ordering_matrix`, not the labelling of the nodes. + tuple[frozenset[int], ...] + Topological generations. Integers represent the indices of the matrix `ordering_matrix`, not the labelling of the nodes. or `None` if `ordering_matrix` is not a DAG. @@ -558,31 +559,16 @@ def _compute_topological_generations(ordering_matrix: MatGF2) -> list[list[int]] adj_mat = ordering_matrix indegree_map: dict[int, int] = {} - zero_indegree: list[int] = [] + zero_indegree: set[int] = set() neighbors = {node: set(np.flatnonzero(row).astype(int)) for node, row in enumerate(adj_mat.T)} for node, col in enumerate(adj_mat): parents = np.flatnonzero(col) if parents.size: indegree_map[node] = parents.size else: - zero_indegree.append(node) - - generations: list[list[int]] = [] - - while zero_indegree: - this_generation = zero_indegree - zero_indegree = [] - for node in this_generation: - for child in neighbors[node]: - indegree_map[child] -= 1 - if indegree_map[child] == 0: - zero_indegree.append(child) - del indegree_map[child] - generations.append(this_generation) - - if indegree_map: - return None - return generations + zero_indegree.add(node) + + return compute_topological_generations(neighbors, indegree_map, zero_indegree) def compute_partial_order_layers(correction_matrix: CorrectionMatrix[_M_co]) -> tuple[frozenset[int], ...] | None: @@ -612,7 +598,7 @@ def compute_partial_order_layers(correction_matrix: CorrectionMatrix[_M_co]) -> aog, c_matrix = correction_matrix.aog, correction_matrix.c_matrix ordering_matrix = aog.order_demand_matrix.mat_mul(c_matrix) - if (topo_gen := _compute_topological_generations(ordering_matrix)) is None: + if (topo_gen := _try_ordering_matrix_to_topological_generations(ordering_matrix)) is None: return None # The NC matrix is not a DAG, therefore there's no flow. layers = [ diff --git a/graphix/flow/_partial_order.py b/graphix/flow/_partial_order.py new file mode 100644 index 000000000..569419c29 --- /dev/null +++ b/graphix/flow/_partial_order.py @@ -0,0 +1,54 @@ +"""Tools for computing the partial orders.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from collections.abc import Mapping + from collections.abc import Set as AbstractSet + + +def compute_topological_generations( + dag: Mapping[int, AbstractSet[int]], indegree_map: Mapping[int, int], zero_indegree: AbstractSet[int] +) -> tuple[frozenset[int], ...] | None: + """Stratify the directed acyclic graph (DAG) into generations. + + Parameters + ---------- + dag : Mapping[int, AbstractSet[int]] + Mapping encoding the directed acyclic graph. + + indegree_map : Mapping[int, int] + Indegree of the input DAG. A pair (``key``, ``value``) represents a node in the DAG and the number of incoming edges incident on it. It is assumed that indegree values are larger than 0. + + zero_indegree : AbstractSet[int] + Nodes in the DAG without any incoming edges. + + Returns + ------- + tuple[frozenset[int], ...] | None + Topological generations. ``None`` if the input DAG contains closed loops. + + Notes + ----- + This function is adapted from `:func: networkx.algorithms.dag.topological_generations` so that it works directly on dictionaries instead of a `:class: nx.DiGraph` object. + """ + generations: list[frozenset[int]] = [] + indegree_map = dict(indegree_map) + + while zero_indegree: + this_generation = zero_indegree + zero_indegree = set() + for node in this_generation: + for child in dag[node]: + indegree_map[child] -= 1 + if indegree_map[child] == 0: + zero_indegree.add(child) + del indegree_map[child] + generations.append(frozenset(this_generation)) + + if indegree_map: + return None + + return tuple(generations) diff --git a/graphix/optimization.py b/graphix/optimization.py index cc67f2833..893debe30 100644 --- a/graphix/optimization.py +++ b/graphix/optimization.py @@ -3,6 +3,7 @@ from __future__ import annotations import dataclasses +from collections import defaultdict from copy import copy from dataclasses import dataclass from types import MappingProxyType @@ -17,8 +18,16 @@ from graphix import command from graphix.clifford import Clifford from graphix.command import CommandKind, Node +from graphix.flow._partial_order import compute_topological_generations +from graphix.flow.core import CausalFlow, GFlow +from graphix.flow.exceptions import ( + FlowGenericError, + FlowGenericErrorReason, +) from graphix.fundamentals import Axis, Plane -from graphix.measurements import Domains, Outcome, PauliMeasurement +from graphix.measurements import Domains, Measurement, Outcome, PauliMeasurement +from graphix.opengraph import OpenGraph +from graphix.states import BasicStates if TYPE_CHECKING: from collections.abc import Iterable, Mapping @@ -362,6 +371,180 @@ def ensure_neighborhood(node: Node) -> None: done.add(node) return pattern + def extract_opengraph(self) -> OpenGraph[Measurement]: + """Extract the underlying resource-state open graph from the pattern. + + Returns + ------- + OpenGraph[Measurement] + + Raises + ------ + ValueError + If `N` commands in the pattern do not represent a |+⟩ state. + + Notes + ----- + This operation loses all the information on the Clifford commands. + """ + for n in self.n_list: + if n.state != BasicStates.PLUS: + raise ValueError( + f"Open graph construction in flow extraction requires N commands to represent a |+⟩ state. Error found in {n}." + ) + measurements = {m.node: Measurement(m.angle, m.plane) for m in self.m_list} + return OpenGraph(self.extract_graph(), self.input_nodes, self.output_nodes, measurements) + + def extract_partial_order_layers(self) -> tuple[frozenset[int], ...]: + """Extract the measurement order of the pattern in the form of layers. + + This method builds a directed acyclical graph (DAG) from the pattern and then performs a topological sort. + + Returns + ------- + tuple[frozenset[int], ...] + Measurement partial order between the pattern's nodes in a layer form. + + Raises + ------ + ValueError + If the correction domains in the pattern form closed loops. This implies that the pattern is not runnable. + + Notes + ----- + The returned object follows the same conventions as the ``partial_order_layers`` attribute of :class:`PauliFlow` and :class:`XZCorrections` objects: + - Nodes in the same layer can be measured simultaneously. + - Nodes in layer ``i`` must be measured after nodes in layer ``i + 1``. + - All output nodes (if any) are in the first layer. + - There cannot be any empty layers. + """ + oset = frozenset(self.output_nodes) # First layer by convention. + pre_measured_nodes = self.results.keys() # Not included in the partial order layers. + excluded_nodes = oset | pre_measured_nodes + + zero_indegree = set(self.input_nodes).union(n.node for n in self.n_list) - excluded_nodes + dag: dict[int, set[int]] = { + node: set() for node in zero_indegree + } # `i: {j}` represents `i -> j` which means that node `i` must be measured before node `j`. + indegree_map: dict[int, int] = {} + + def process_domain(node: Node, domain: AbstractSet[Node]) -> None: + for dep_node in domain: + if ( + not {node, dep_node} & excluded_nodes and node not in dag[dep_node] + ): # Don't include multiple edges in the dag. + dag[dep_node].add(node) + indegree_map[node] = indegree_map.get(node, 0) + 1 + + domain: AbstractSet[Node] + + for m in self.m_list: + node, domain = m.node, m.s_domain | m.t_domain + process_domain(node, domain) + + for corrections in [self.z_dict, self.x_dict]: + for node, domain in corrections.items(): + process_domain(node, domain) + + zero_indegree -= indegree_map.keys() + + generations = compute_topological_generations(dag, indegree_map, zero_indegree) + if generations is None: + raise ValueError("Pattern domains form closed loops.") + + if oset: + return oset, *generations[::-1] + return generations[::-1] + + def extract_causal_flow(self) -> CausalFlow[Measurement]: + """Extract the causal flow structure from the current measurement pattern. + + This method does not call the flow-extraction routine on the underlying open graph, but constructs the flow from the pattern corrections instead. + + Returns + ------- + flow.CausalFlow[Measurement] + The causal flow associated with the current pattern. + + Raises + ------ + FlowError + If the pattern: + - Contains measurements in forbidden planes (XZ or YZ), + - Is empty, or + - Induces a correction function and a partial order which fail the well-formedness checks for a valid causal flow. + ValueError + If `N` commands in the pattern do not represent a |+⟩ state or if the pattern corrections form closed loops. + + Notes + ----- + This method makes use of :func:`StandardizedPattern.extract_partial_order_layers` which computes the pattern's direct acyclical graph (DAG) induced by the corrections and returns a particular layer stratification (obtained by doing a topological sort on the DAG). Further, it constructs the pattern's induced correction function from :math:`M` and :math:`X` commands. + In general, there may exist various layerings which represent the corrections of the pattern. To ensure that a given layering is compatible with the pattern's induced correction function, the partial order must be extracted from a standardized pattern. Commutation of entanglement commands with X and Z corrections in the standardization procedure may generate new corrections, which guarantees that all the topological information of the underlying graph is encoded in the extracted partial order. + """ + correction_function: dict[int, set[int]] = defaultdict(set) + pre_measured_nodes = self.results.keys() # Not included in the flow. + + for m in self.m_list: + if m.plane in {Plane.XZ, Plane.YZ}: + raise FlowGenericError(FlowGenericErrorReason.XYPlane) + _update_corrections(m.node, m.s_domain - pre_measured_nodes, correction_function) + + for node, domain in self.x_dict.items(): + _update_corrections(node, domain - pre_measured_nodes, correction_function) + + og = ( + self.extract_opengraph() + ) # Raises a `ValueError` if `N` commands in the pattern do not represent a |+⟩ state. + partial_order_layers = ( + self.extract_partial_order_layers() + ) # Raises a `ValueError` if the pattern corrections form closed loops. + cf = CausalFlow(og, dict(correction_function), partial_order_layers) + cf.check_well_formed() # Raises a `FlowError` if the partial order and the correction function are not compatible, or if a measured node is corrected by more than one node. + return cf + + def extract_gflow(self) -> GFlow[Measurement]: + """Extract the generalized flow (gflow) structure from the current measurement pattern. + + This method does not call the flow-extraction routine on the underlying open graph, but constructs the gflow from the pattern corrections instead. + + Returns + ------- + flow.GFlow[Measurement] + The gflow associated with the current pattern. + + Raises + ------ + FlowError + If the pattern is empty or if the extracted structure does not satisfy + the well-formedness conditions required for a valid gflow. + ValueError + If `N` commands in the pattern do not represent a |+⟩ state or if the pattern corrections form closed loops. + + Notes + ----- + The notes provided in :func:`self.extract_causal_flow` apply here as well. + """ + correction_function: dict[int, set[int]] = {} + pre_measured_nodes = self.results.keys() # Not included in the flow. + + for m in self.m_list: + if m.plane in {Plane.XZ, Plane.YZ}: + correction_function.setdefault(m.node, set()).add(m.node) + _update_corrections(m.node, m.s_domain - pre_measured_nodes, correction_function) + + for node, domain in self.x_dict.items(): + _update_corrections(node, domain - pre_measured_nodes, correction_function) + + og = ( + self.extract_opengraph() + ) # Raises a `ValueError` if `N` commands in the pattern do not represent a |+⟩ state. + partial_order_layers = ( + self.extract_partial_order_layers() + ) # Raises a `ValueError` if the pattern corrections form closed loops. + gf = GFlow(og, correction_function, partial_order_layers) + gf.check_well_formed() + return gf + def _add_correction_domain(domain_dict: dict[Node, set[Node]], node: Node, domain: set[Node]) -> None: """Merge a correction domain into ``domain_dict`` for ``node``. @@ -418,6 +601,27 @@ def _incorporate_pauli_results_in_domain( return odd_outcome == 1, new_domain +def _update_corrections(node: Node, domain: AbstractSet[Node], correction: dict[Node, set[Node]]) -> None: + """Update the correction mapping by adding a node to all entries in a domain. + + Parameters + ---------- + node : Node + The node to add as a correction. + domain : AbstractSet[Node] + A set of measured nodes whose corresponding correction sets should be updated. + correction : dict[Node, set[Node]] + A mapping from measured nodes to sets of nodes on which corrections are applied. This + dictionary is modified in place. + + Notes + ----- + This function is used to extract the correction function from :math:`X`, :math:`Z` and :math:`M` commands when constructing a flow. + """ + for measured_node in domain: + correction.setdefault(measured_node, set()).add(node) + + def incorporate_pauli_results(pattern: Pattern) -> Pattern: """Return an equivalent pattern where results from Pauli presimulation are integrated in corrections.""" result = graphix.pattern.Pattern(input_nodes=pattern.input_nodes) diff --git a/graphix/pattern.py b/graphix/pattern.py index dbc312553..4ff84c602 100644 --- a/graphix/pattern.py +++ b/graphix/pattern.py @@ -8,6 +8,7 @@ import copy import dataclasses import enum +import itertools import warnings from collections.abc import Iterable, Iterator from dataclasses import dataclass @@ -22,7 +23,6 @@ from graphix.clifford import Clifford from graphix.command import Command, CommandKind from graphix.fundamentals import Axis, Plane, Sign -from graphix.gflow import find_flow, find_gflow, get_layers from graphix.graphsim import GraphState from graphix.measurements import Measurement, Outcome, PauliMeasurement, toggle_outcome from graphix.opengraph import OpenGraph @@ -38,6 +38,7 @@ from numpy.random import Generator + from graphix.flow.core import CausalFlow, GFlow from graphix.parameter import ExpressionOrFloat, ExpressionOrSupportsFloat, Parameter from graphix.sim import Backend, BackendState, Data @@ -885,6 +886,80 @@ def update_dependency(measured: AbstractSet[int], dependency: dict[int, set[int] for i in dependency: dependency[i] -= measured + def extract_partial_order_layers(self) -> tuple[frozenset[int], ...]: + """Extract the measurement order of the pattern in the form of layers. + + This method standardizes the pattern, builds a directed acyclical graph (DAG) from measurement and correction domains, and then performs a topological sort. + + Returns + ------- + tuple[frozenset[int], ...] + Measurement partial order between the pattern's nodes in a layer form. + + Raises + ------ + RunnabilityError + If the pattern is not runnable. + + Notes + ----- + - This function wraps :func:`optimization.StandardizedPattern.extract_partial_order_layers`, and the returned object is described in the notes of this method. + + - See :func:`optimization.StandardizedPattern.extract_causal_flow` for additional information on why it is required to standardized the pattern to extract the partial order layering. + """ + return optimization.StandardizedPattern.from_pattern(self).extract_partial_order_layers() + + def extract_causal_flow(self) -> CausalFlow[Measurement]: + """Extract the causal flow structure from the current measurement pattern. + + This method does not call the flow-extraction routine on the underlying open graph, but constructs the flow from the pattern corrections instead. + + Returns + ------- + CausalFlow[Measurement] + The causal flow associated with the current pattern. + + Raises + ------ + FlowError + If the pattern: + - Contains measurements in forbidden planes (XZ or YZ), + - Is empty, or + - Induces a correction function and a partial order which fail the well-formedness checks for a valid causal flow. + ValueError + If `N` commands in the pattern do not represent a |+⟩ state or if the pattern corrections form closed loops. + + Notes + ----- + - See :func:`optimization.StandardizedPattern.extract_causal_flow` for additional information on why it is required to standardized the pattern to extract a causal flow. + - Applying the chain ``Pattern.extract_causal_flow().to_corrections().to_pattern()`` to a strongly deterministic pattern returns a new pattern implementing the same unitary transformation. This equivalence holds as long as the original pattern contains no Clifford commands, since those are discarded during open-graph extraction. + """ + return optimization.StandardizedPattern.from_pattern(self).extract_causal_flow() + + def extract_gflow(self) -> GFlow[Measurement]: + """Extract the generalized flow (gflow) structure from the current measurement pattern. + + This method does not call the flow-extraction routine on the underlying open graph, but constructs the gflow from the pattern corrections instead. + + Returns + ------- + GFlow[Measurement] + The gflow associated with the current pattern. + + Raises + ------ + FlowError + If the pattern is empty or if the extracted structure does not satisfy + the well-formedness conditions required for a valid gflow. + ValueError + If `N` commands in the pattern do not represent a |+⟩ state or if the pattern corrections form closed loops. + + Notes + ----- + The notes provided in :func:`self.extract_causal_flow` apply here as well. + """ + return optimization.StandardizedPattern.from_pattern(self).extract_gflow() + def get_layers(self) -> tuple[int, dict[int, set[int]]]: """Construct layers(l_k) from dependency information. @@ -925,14 +1000,11 @@ def _measurement_order_depth(self) -> list[int]: Returns ------- - meas_order: list of int + list[int] optimal measurement order for parallel computing """ - d, l_k = self.get_layers() - meas_order: list[int] = [] - for i in range(d): - meas_order.extend(l_k[i]) - return meas_order + partial_order_layers = self.extract_partial_order_layers() + return list(itertools.chain(*reversed(partial_order_layers[1:]))) @staticmethod def connected_edges(node: int, edges: set[tuple[int, int]]) -> set[tuple[int, int]]: @@ -984,54 +1056,6 @@ def _measurement_order_space(self) -> list[int]: edges -= removable_edges return meas_order - def get_measurement_order_from_flow(self) -> list[int] | None: - """Return a measurement order generated from flow. If a graph has flow, the minimum 'max_space' of a pattern is guaranteed to width+1. - - Returns - ------- - meas_order: list of int - measurement order - """ - graph = self.extract_graph() - vin = set(self.input_nodes) - vout = set(self.output_nodes) - meas_planes = self.get_meas_plane() - f, l_k = find_flow(graph, vin, vout, meas_planes=meas_planes) - if f is None or l_k is None: - return None - depth, layer = get_layers(l_k) - meas_order: list[int] = [] - for i in range(depth): - k = depth - i - nodes = layer[k] - meas_order += nodes # NOTE this is list concatenation - return meas_order - - def get_measurement_order_from_gflow(self) -> list[int]: - """Return a list containing the node indices, in the order of measurements which can be performed with minimum depth. - - Returns - ------- - meas_order : list of int - measurement order - """ - graph = self.extract_graph() - isolated = list(nx.isolates(graph)) - if isolated: - raise ValueError("The input graph must be connected") - vin = set(self.input_nodes) - vout = set(self.output_nodes) - meas_planes = self.get_meas_plane() - flow, l_k = find_gflow(graph, vin, vout, meas_planes=meas_planes) - if flow is None or l_k is None: # We check both to avoid typing issues with `get_layers`. - raise ValueError("No gflow found") - k, layers = get_layers(l_k) - meas_order: list[int] = [] - while k > 0: - meas_order.extend(layers[k]) - k -= 1 - return meas_order - def sort_measurement_commands(self, meas_order: list[int]) -> list[command.M]: """Convert measurement order to sequence of measurement commands. @@ -1287,7 +1311,8 @@ def minimize_space(self) -> None: self.standardize() meas_order = None if not self._pauli_preprocessed: - meas_order = self.get_measurement_order_from_flow() + cf = self.extract_opengraph().find_causal_flow() + meas_order = list(itertools.chain(*reversed(cf.partial_order_layers[1:]))) if cf is not None else None if meas_order is None: meas_order = self._measurement_order_space() self._reorder_pattern(self.sort_measurement_commands(meas_order)) diff --git a/graphix/visualization.py b/graphix/visualization.py index c2e938707..ffcd295b6 100644 --- a/graphix/visualization.py +++ b/graphix/visualization.py @@ -11,8 +11,11 @@ from matplotlib import pyplot as plt from graphix import gflow +from graphix.flow.exceptions import FlowError from graphix.fundamentals import Plane from graphix.measurements import PauliMeasurement +from graphix.opengraph import OpenGraph +from graphix.optimization import StandardizedPattern if TYPE_CHECKING: from collections.abc import Callable, Collection, Hashable, Iterable, Mapping, Sequence @@ -24,6 +27,8 @@ # MEMO: Potential circular import from graphix.clifford import Clifford + from graphix.flow.core import CausalFlow, GFlow + from graphix.fundamentals import AbstractPlanarMeasurement from graphix.parameter import ExpressionOrFloat from graphix.pattern import Pattern @@ -127,28 +132,38 @@ def visualize( If not None, filename of the png file to save the plot. If None, the plot is not saved. Default in None. """ - f, l_k = gflow.find_flow(self.graph, set(self.v_in), set(self.v_out), meas_planes=self.meas_planes) # try flow - if f is not None and l_k is not None: + og = OpenGraph(self.graph, list(self.v_in), list(self.v_out), self.meas_planes) + causal_flow = og.find_causal_flow() + if causal_flow is not None: print("Flow detected in the graph.") - pos = self.get_pos_from_flow(f, l_k) + pos = self.get_pos_from_flow(causal_flow) + cf = causal_flow.correction_function + l_k = { + node: layer_idx for layer_idx, layer in enumerate(causal_flow.partial_order_layers) for node in layer + } def get_paths( pos: Mapping[int, _Point], ) -> tuple[Mapping[_Edge, Sequence[_Point]], Mapping[_Edge, Sequence[_Point]] | None]: - return self.get_edge_path(f, pos) + return self.get_edge_path(cf, pos) + else: - g, l_k = gflow.find_gflow(self.graph, set(self.v_in), set(self.v_out), self.meas_planes) # try gflow - if g is not None and l_k is not None: + g_flow = og.find_gflow() + if g_flow is not None: print("Gflow detected in the graph. (flow not detected)") - pos = self.get_pos_from_gflow(g, l_k) + pos = self.get_pos_from_gflow(g_flow) + cf = g_flow.correction_function + l_k = {node: layer_idx for layer_idx, layer in enumerate(g_flow.partial_order_layers) for node in layer} def get_paths( pos: Mapping[int, _Point], ) -> tuple[Mapping[_Edge, Sequence[_Point]], Mapping[_Edge, Sequence[_Point]] | None]: - return self.get_edge_path(g, pos) + return self.get_edge_path(cf, pos) + else: print("No flow or gflow detected in the graph.") pos = self.get_pos_wo_structure() + l_k = None def get_paths( pos: Mapping[int, _Point], @@ -207,35 +222,19 @@ def visualize_from_pattern( If not None, filename of the png file to save the plot. If None, the plot is not saved. Default in None. """ - f, l_k = gflow.flow_from_pattern(pattern) # try flow - if f is not None and l_k is not None: - print("The pattern is consistent with flow structure.") - pos = self.get_pos_from_flow(f, l_k) - - def get_paths( - pos: Mapping[int, _Point], - ) -> tuple[Mapping[_Edge, Sequence[_Point]], Mapping[_Edge, Sequence[_Point]] | None]: - return self.get_edge_path(f, pos) - - corrections: tuple[Mapping[int, AbstractSet[int]], Mapping[int, AbstractSet[int]]] | None = None - else: - g, l_k = gflow.gflow_from_pattern(pattern) # try gflow - if g is not None and l_k is not None: - print("The pattern is consistent with gflow structure. (not with flow)") - pos = self.get_pos_from_gflow(g, l_k) - - def get_paths( - pos: Mapping[int, _Point], - ) -> tuple[Mapping[_Edge, Sequence[_Point]], Mapping[_Edge, Sequence[_Point]] | None]: - return self.get_edge_path(g, pos) + pattern_std = StandardizedPattern.from_pattern(pattern) + cf: Mapping[int, AbstractSet[int]] + corrections: tuple[Mapping[int, AbstractSet[int]], Mapping[int, AbstractSet[int]]] | None - corrections = None - else: + try: + causal_flow = pattern_std.extract_causal_flow() + except FlowError: + try: + g_flow = pattern_std.extract_gflow() + except FlowError: print("The pattern is not consistent with flow or gflow structure.") - depth, layers = pattern.get_layers() - unfolded_layers = {element: key for key, value_set in layers.items() for element in value_set} - for output in pattern.output_nodes: - unfolded_layers[output] = depth + 1 + po_layers = pattern.extract_partial_order_layers() + unfolded_layers = {node: layer_idx for layer_idx, layer in enumerate(po_layers[::-1]) for node in layer} xflow, zflow = gflow.get_corrections_from_pattern(pattern) xzflow: dict[int, set[int]] = deepcopy(xflow) for key, value in zflow.items(): @@ -244,13 +243,29 @@ def get_paths( else: xzflow[key] = set(value) # copy pos = self.get_pos_all_correction(unfolded_layers) - - def get_paths( - pos: Mapping[int, _Point], - ) -> tuple[Mapping[_Edge, Sequence[_Point]], Mapping[_Edge, Sequence[_Point]] | None]: - return self.get_edge_path(xzflow, pos) - + cf = xzflow + l_k = None corrections = xflow, zflow + else: + print("The pattern is consistent with gflow structure. (not with flow)") + pos = self.get_pos_from_gflow(g_flow) + cf = g_flow.correction_function + l_k = {node: layer_idx for layer_idx, layer in enumerate(g_flow.partial_order_layers) for node in layer} + corrections = None + else: + print("The pattern is consistent with flow structure.") + pos = self.get_pos_from_flow(causal_flow) + cf = causal_flow.correction_function + l_k = { + node: layer_idx for layer_idx, layer in enumerate(causal_flow.partial_order_layers) for node in layer + } + corrections = None + + def get_paths( + pos: Mapping[int, _Point], + ) -> tuple[Mapping[_Edge, Sequence[_Point]], Mapping[_Edge, Sequence[_Point]] | None]: + return self.get_edge_path(cf, pos) + self.visualize_graph( pos, get_paths, @@ -511,7 +526,7 @@ def get_figsize( return (width * node_distance[0], height * node_distance[1]) def get_edge_path( - self, flow: Mapping[int, int | set[int]], pos: Mapping[int, _Point] + self, flow: Mapping[int, AbstractSet[int]], pos: Mapping[int, _Point] ) -> tuple[dict[_Edge, list[_Point]], dict[_Edge, list[_Point]]]: """ Return the path of edges and gflow arrows. @@ -533,7 +548,7 @@ def get_edge_path( edge_path = self.get_edge_path_wo_structure(pos) edge_set = set(self.graph.edges()) arrow_path: dict[_Edge, list[_Point]] = {} - flow_arrows = {(k, v) for k, values in flow.items() for v in ((values,) if isinstance(values, int) else values)} + flow_arrows = {(k, v) for k, values in flow.items() for v in values} for arrow in flow_arrows: if arrow[0] == arrow[1]: # Self loop @@ -632,7 +647,7 @@ def get_edge_path_wo_structure(self, pos: Mapping[int, _Point]) -> dict[_Edge, l """ return {edge: self._find_bezier_path(edge, [pos[edge[0]], pos[edge[1]]], pos) for edge in self.graph.edges()} - def get_pos_from_flow(self, f: Mapping[int, set[int]], l_k: Mapping[int, int]) -> dict[int, _Point]: + def get_pos_from_flow(self, flow: CausalFlow[AbstractPlanarMeasurement]) -> dict[int, _Point]: """ Return the position of nodes based on the flow. @@ -648,6 +663,7 @@ def get_pos_from_flow(self, f: Mapping[int, set[int]], l_k: Mapping[int, int]) - pos : dict dictionary of node positions. """ + f = flow.correction_function values_union = set().union(*f.values()) start_nodes = set(self.graph.nodes()) - values_union pos = {node: [0, 0] for node in self.graph.nodes()} @@ -658,16 +674,15 @@ def get_pos_from_flow(self, f: Mapping[int, set[int]], l_k: Mapping[int, int]) - node = next(iter(f[node])) pos[node][1] = i - if not l_k: - return {} - - lmax = max(l_k.values()) + layers = flow.partial_order_layers + lmax = len(layers) - 1 # Change the x coordinates of the nodes based on their layer, sort in descending order - for node, layer in l_k.items(): - pos[node][0] = lmax - layer + for layer_idx, layer in enumerate(layers): + for node in layer: + pos[node][0] = lmax - layer_idx return {k: (x, y) for k, (x, y) in pos.items()} - def get_pos_from_gflow(self, g: Mapping[int, set[int]], l_k: Mapping[int, int]) -> dict[int, _Point]: + def get_pos_from_gflow(self, flow: GFlow[AbstractPlanarMeasurement]) -> dict[int, _Point]: """ Return the position of nodes based on the gflow. @@ -685,6 +700,8 @@ def get_pos_from_gflow(self, g: Mapping[int, set[int]], l_k: Mapping[int, int]) """ g_edges: list[_Edge] = [] + g = flow.correction_function + for node, node_list in g.items(): g_edges.extend((node, n) for n in node_list) @@ -692,17 +709,19 @@ def get_pos_from_gflow(self, g: Mapping[int, set[int]], l_k: Mapping[int, int]) g_prime.add_nodes_from(self.graph.nodes()) g_prime.add_edges_from(g_edges) - l_max = max(l_k.values()) - l_reverse = {v: l_max - l for v, l in l_k.items()} + layers = flow.partial_order_layers + l_max = len(layers) - 1 + l_reverse = {node: l_max - layer_idx for layer_idx, layer in enumerate(layers) for node in layer} _set_node_attributes(g_prime, l_reverse, "subset") - pos = nx.multipartite_layout(g_prime) vert = list({pos[node][1] for node in self.graph.nodes()}) vert.sort() index = {y: i for i, y in enumerate(vert)} - return {node: (l_max - layer, index[pos[node][1]]) for node, layer in l_k.items()} + return { + node: (l_max - layer_idx, index[pos[node][1]]) for layer_idx, layer in enumerate(layers) for node in layer + } def get_pos_wo_structure(self) -> dict[int, _Point]: """ diff --git a/tests/baseline/test_draw_graph_reference_True.png b/tests/baseline/test_draw_graph_reference_True.png index 27f65ee5a..2f79065a6 100644 Binary files a/tests/baseline/test_draw_graph_reference_True.png and b/tests/baseline/test_draw_graph_reference_True.png differ diff --git a/tests/test_flow_find_gpflow.py b/tests/test_flow_find_gpflow.py index 704008db2..c0dfeda25 100644 --- a/tests/test_flow_find_gpflow.py +++ b/tests/test_flow_find_gpflow.py @@ -21,7 +21,7 @@ from graphix.flow._find_gpflow import ( AlgebraicOpenGraph, PlanarAlgebraicOpenGraph, - _compute_topological_generations, + _try_ordering_matrix_to_topological_generations, compute_correction_matrix, ) from graphix.fundamentals import Axis, Plane @@ -42,7 +42,7 @@ class AlgebraicOpenGraphTestCase(NamedTuple): class DAGTestCase(NamedTuple): adj_mat: MatGF2 - generations: list[list[int]] | None + generations: tuple[frozenset[int], ...] | None def prepare_test_og() -> list[AlgebraicOpenGraphTestCase]: @@ -230,14 +230,15 @@ def prepare_test_dag() -> list[DAGTestCase]: test_cases.extend( ( # Simple DAG DAGTestCase( - adj_mat=MatGF2([[0, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0], [0, 1, 1, 0]]), generations=[[0], [1, 2], [3]] + adj_mat=MatGF2([[0, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0], [0, 1, 1, 0]]), + generations=(frozenset({0}), frozenset({1, 2}), frozenset({3})), ), # Graph with loop DAGTestCase(adj_mat=MatGF2([[0, 0, 0, 0], [1, 0, 0, 1], [1, 0, 0, 0], [0, 1, 1, 0]]), generations=None), # Disconnected graph DAGTestCase( adj_mat=MatGF2([[0, 0, 0, 0, 0], [1, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 1, 0, 0]]), - generations=[[0, 2], [1, 3, 4]], + generations=(frozenset({0, 2}), frozenset({1, 3, 4})), ), ) ) @@ -273,8 +274,8 @@ def test_correction_matrix(self, test_case: AlgebraicOpenGraphTestCase) -> None: assert corr_matrix is None @pytest.mark.parametrize("test_case", prepare_test_dag()) - def test_compute_topological_generations(self, test_case: DAGTestCase) -> None: + def test_try_ordering_matrix_to_topological_generations(self, test_case: DAGTestCase) -> None: adj_mat = test_case.adj_mat generations_ref = test_case.generations - assert generations_ref == _compute_topological_generations(adj_mat) + assert generations_ref == _try_ordering_matrix_to_topological_generations(adj_mat) diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 7693845df..1e2fdf18b 100644 --- a/tests/test_optimization.py +++ b/tests/test_optimization.py @@ -7,7 +7,6 @@ from graphix.clifford import Clifford from graphix.command import C, Command, CommandKind, E, M, N, X, Z from graphix.fundamentals import Plane -from graphix.gflow import gflow_from_pattern from graphix.optimization import StandardizedPattern, incorporate_pauli_results, remove_useless_domains from graphix.pattern import Pattern from graphix.random_objects import rand_circuit @@ -81,9 +80,8 @@ def test_flow_after_pauli_preprocessing(fx_bg: PCG64, jumps: int) -> None: # pattern.move_pauli_measurements_to_the_front() pattern.perform_pauli_measurements() pattern2 = incorporate_pauli_results(pattern) - pattern2.standardize() - f, _l = gflow_from_pattern(pattern2) - assert f is not None + gflow = pattern2.extract_gflow() + gflow.check_well_formed() @pytest.mark.parametrize("jumps", range(1, 11)) diff --git a/tests/test_pattern.py b/tests/test_pattern.py index fc4a61b5a..eadbcb6e6 100644 --- a/tests/test_pattern.py +++ b/tests/test_pattern.py @@ -3,7 +3,7 @@ import copy import itertools import typing -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, NamedTuple import networkx as nx import numpy as np @@ -13,8 +13,12 @@ from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector from graphix.clifford import Clifford from graphix.command import C, Command, CommandKind, E, M, N, X, Z +from graphix.flow.exceptions import ( + FlowError, +) from graphix.fundamentals import Plane -from graphix.measurements import Outcome, PauliMeasurement +from graphix.measurements import Measurement, Outcome, PauliMeasurement +from graphix.opengraph import OpenGraph from graphix.pattern import Pattern, RunnabilityError, RunnabilityErrorReason, shift_outcomes from graphix.random_objects import rand_circuit, rand_gate from graphix.sim.density_matrix import DensityMatrix @@ -751,7 +755,7 @@ def test_check_runnability_failures(self) -> None: pattern = Pattern(cmds=[N(0), M(0, s_domain={0})]) with pytest.raises(RunnabilityError) as exc_info: - pattern.get_layers() + pattern.extract_partial_order_layers() assert exc_info.value.node == 0 assert exc_info.value.reason == RunnabilityErrorReason.DomainSelfLoop @@ -770,6 +774,260 @@ def test_check_runnability_failures(self) -> None: def test_compute_max_degree_empty_pattern(self) -> None: assert Pattern().compute_max_degree() == 0 + @pytest.mark.parametrize( + "test_case", + [ + ( + Pattern(input_nodes=[0], cmds=[N(1), E((0, 1)), M(0), M(1)]), + (frozenset({0, 1}),), + ), + ( + Pattern(input_nodes=[0], cmds=[N(1), N(2), E((0, 1)), E((1, 2)), M(0), M(1), X(2, {1}), Z(2, {0})]), + (frozenset({2}), frozenset({0, 1})), + ), + ( + Pattern(input_nodes=[0, 1], cmds=[M(1), M(0, s_domain={1}), N(2)]), + (frozenset({2}), frozenset({0}), frozenset({1})), + ), + ( + Pattern( + input_nodes=[0], cmds=[N(1), N(2), E((0, 1)), E((1, 2)), M(0), M(1), X(2, {1}), Z(2, {1}), M(2)] + ), + (frozenset({2}), frozenset({0, 1})), + ), # double edge in DAG + ], + ) + def test_extract_partial_order_layers(self, test_case: tuple[Pattern, tuple[frozenset[int], ...]]) -> None: + assert test_case[0].extract_partial_order_layers() == test_case[1] + + def test_extract_partial_order_layers_results(self) -> None: + c = Circuit(1) + c.rz(0, 0.2) + p = c.transpile().pattern + p.perform_pauli_measurements() + assert p.extract_partial_order_layers() == (frozenset({2}), frozenset({0})) + + p = Pattern(cmds=[N(0), N(1), N(2), M(0), E((1, 2)), X(1, {0}), M(2, angle=0.3)]) + p.perform_pauli_measurements() + assert p.extract_partial_order_layers() == (frozenset({1}), frozenset({2})) + + class PatternFlowTestCase(NamedTuple): + pattern: Pattern + has_cflow: bool + has_gflow: bool + + PATTERN_FLOW_TEST_CASES: list[PatternFlowTestCase] = [ # noqa: RUF012 + PatternFlowTestCase( + # General example + Pattern( + input_nodes=[0, 1], + cmds=[ + N(2), + N(3), + N(4), + N(5), + N(6), + N(7), + E((0, 2)), + E((2, 3)), + E((2, 4)), + E((1, 3)), + E((3, 5)), + E((4, 5)), + E((4, 6)), + E((5, 7)), + M(0, angle=0.1), + Z(3, {0}), + Z(4, {0}), + X(2, {0}), + M(1, angle=0.1), + Z(2, {1}), + Z(5, {1}), + X(3, {1}), + M(2, angle=0.1), + Z(5, {2}), + Z(6, {2}), + X(4, {2}), + M(3, angle=0.1), + Z(4, {3}), + Z(7, {3}), + X(5, {3}), + M(4, angle=0.1), + X(6, {4}), + M(5, angle=0.4), + X(7, {5}), + ], + output_nodes=[6, 7], + ), + has_cflow=True, + has_gflow=True, + ), + PatternFlowTestCase( + # No measurements or corrections + Pattern(input_nodes=[0, 1], cmds=[E((0, 1))]), + has_cflow=True, + has_gflow=True, + ), + PatternFlowTestCase( + # Disconnected nodes and unordered outputs + Pattern(input_nodes=[2], cmds=[N(0), N(1), E((0, 1)), M(0), X(1, {0})], output_nodes=[2, 1]), + has_cflow=True, + has_gflow=True, + ), + PatternFlowTestCase( + # Pattern with XZ measurements. + Pattern(cmds=[N(0), N(1), E((0, 1)), M(0, Plane.XZ, 0.3), Z(1, {0}), X(1, {0})], output_nodes=[1]), + has_cflow=False, + has_gflow=True, + ), + PatternFlowTestCase( + # Pattern with gflow but without causal flow and XY measurements. + Pattern( + input_nodes=[1, 2, 3], + cmds=[ + N(4), + N(5), + N(6), + E((1, 4)), + E((1, 6)), + E((4, 2)), + E((6, 2)), + E((6, 3)), + E((2, 5)), + E((5, 3)), + M(1, angle=0.1), + X(5, {1}), + X(6, {1}), + M(2, angle=0.2), + X(4, {2}), + X(5, {2}), + X(6, {2}), + M(3, angle=0.3), + X(4, {3}), + X(6, {3}), + ], + output_nodes=[4, 5, 6], + ), + has_cflow=False, + has_gflow=True, + ), + PatternFlowTestCase( + # Non-deterministic pattern + Pattern(input_nodes=[0], cmds=[N(1), E((0, 1)), M(0, Plane.XY, 0.3)]), + has_cflow=False, + has_gflow=False, + ), + ] + + # Extract causal flow from random circuits + @pytest.mark.parametrize("jumps", range(1, 11)) + def test_extract_causal_flow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + nqubits = 2 + depth = 2 + circuit_1 = rand_circuit(nqubits, depth, rng, use_ccx=False) + p_ref = circuit_1.transpile().pattern + p_test = p_ref.extract_causal_flow().to_corrections().to_pattern() + + p_ref.perform_pauli_measurements() + p_test.perform_pauli_measurements() + + s_ref = p_ref.simulate_pattern(rng=rng) + s_test = p_test.simulate_pattern(rng=rng) + assert np.abs(np.dot(s_ref.flatten().conjugate(), s_test.flatten())) == pytest.approx(1) + + # Extract gflow from random circuits + @pytest.mark.parametrize("jumps", range(1, 11)) + def test_extract_gflow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None: + rng = Generator(fx_bg.jumped(jumps)) + nqubits = 2 + depth = 2 + circuit_1 = rand_circuit(nqubits, depth, rng, use_ccx=False) + p_ref = circuit_1.transpile().pattern + p_test = p_ref.extract_gflow().to_corrections().to_pattern() + + p_ref.perform_pauli_measurements() + p_test.perform_pauli_measurements() + + s_ref = p_ref.simulate_pattern(rng=rng) + s_test = p_test.simulate_pattern(rng=rng) + assert np.abs(np.dot(s_ref.flatten().conjugate(), s_test.flatten())) == pytest.approx(1) + + @pytest.mark.parametrize("test_case", PATTERN_FLOW_TEST_CASES) + def test_extract_causal_flow(self, fx_rng: Generator, test_case: PatternFlowTestCase) -> None: + if test_case.has_cflow: + alpha = 2 * np.pi * fx_rng.random() + s_ref = test_case.pattern.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha)) + + p_test = test_case.pattern.extract_causal_flow().to_corrections().to_pattern() + s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng) + + assert np.abs(np.dot(s_ref.flatten().conjugate(), s_test.flatten())) == pytest.approx(1) + else: + with pytest.raises(FlowError): + test_case.pattern.extract_causal_flow() + + @pytest.mark.parametrize("test_case", PATTERN_FLOW_TEST_CASES) + def test_extract_gflow(self, fx_rng: Generator, test_case: PatternFlowTestCase) -> None: + if test_case.has_gflow: + alpha = 2 * np.pi * fx_rng.random() + s_ref = test_case.pattern.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha)) + + p_test = test_case.pattern.extract_gflow().to_corrections().to_pattern() + s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng) + + assert np.abs(np.dot(s_ref.flatten().conjugate(), s_test.flatten())) == pytest.approx(1) + else: + with pytest.raises(FlowError): + test_case.pattern.extract_gflow() + + # From open graph + def test_extract_cflow_og(self, fx_rng: Generator) -> None: + alpha = 2 * np.pi * fx_rng.random() + + og = OpenGraph( + graph=nx.Graph([(1, 3), (2, 4), (3, 4), (3, 5), (4, 6)]), + input_nodes=[1, 2], + output_nodes=[6, 5], + measurements={ + 1: Measurement(0.1, Plane.XY), + 2: Measurement(0.2, Plane.XY), + 3: Measurement(0.3, Plane.XY), + 4: Measurement(0.4, Plane.XY), + }, + ) + p_ref = og.extract_causal_flow().to_corrections().to_pattern() + s_ref = p_ref.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha)) + + p_test = p_ref.extract_causal_flow().to_corrections().to_pattern() + s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha)) + + assert np.abs(np.dot(s_ref.flatten().conjugate(), s_test.flatten())) == pytest.approx(1) + + # From open graph + def test_extract_gflow_og(self, fx_rng: Generator) -> None: + alpha = 2 * np.pi * fx_rng.random() + + og = OpenGraph( + graph=nx.Graph([(1, 3), (2, 4), (3, 4), (3, 5), (4, 6)]), + input_nodes=[1, 2], + output_nodes=[6, 5], + measurements={ + 1: Measurement(0.1, Plane.XY), + 2: Measurement(0.2, Plane.XY), + 3: Measurement(0.3, Plane.XY), + 4: Measurement(0.4, Plane.XY), + }, + ) + + p_ref = og.extract_gflow().to_corrections().to_pattern() + s_ref = p_ref.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha)) + + p_test = p_ref.extract_gflow().to_corrections().to_pattern() + s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha)) + + assert np.abs(np.dot(s_ref.flatten().conjugate(), s_test.flatten())) == pytest.approx(1) + def cp(circuit: Circuit, theta: float, control: int, target: int) -> None: """Controlled rotation gate, decomposed.""" # noqa: D401 diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py index dee2fe064..447ed49e9 100644 --- a/tests/test_transpiler.py +++ b/tests/test_transpiler.py @@ -9,7 +9,6 @@ from graphix import instruction from graphix.branch_selector import ConstBranchSelector from graphix.fundamentals import Axis, Sign -from graphix.gflow import flow_from_pattern from graphix.instruction import InstructionKind from graphix.random_objects import rand_circuit, rand_gate, rand_state_vector from graphix.states import BasicStates @@ -235,9 +234,8 @@ def test_add_extend(self) -> None: def test_instruction_flow(self, fx_rng: Generator, instruction: InstructionTestCase) -> None: circuit = Circuit(3, instr=[instruction(fx_rng)]) pattern = circuit.transpile().pattern - pattern.standardize() - f, _l = flow_from_pattern(pattern) - assert f is not None + flow = pattern.extract_causal_flow() + flow.check_well_formed() @pytest.mark.parametrize("jumps", range(1, 11)) @pytest.mark.parametrize("instruction", INSTRUCTION_TEST_CASES) diff --git a/tests/test_visualization.py b/tests/test_visualization.py index e0a9e2120..bc41e2633 100644 --- a/tests/test_visualization.py +++ b/tests/test_visualization.py @@ -9,7 +9,7 @@ import networkx as nx import pytest -from graphix import Circuit, Pattern, command, gflow, visualization +from graphix import Circuit, Pattern, command, visualization from graphix.fundamentals import Plane from graphix.measurements import Measurement from graphix.opengraph import OpenGraph, OpenGraphError @@ -102,10 +102,9 @@ def test_get_pos_from_flow() -> None: meas_angles = pattern.get_angles() local_clifford = pattern.get_vops() vis = visualization.GraphVisualizer(graph, vin, vout, meas_planes, meas_angles, local_clifford) - f, l_k = gflow.find_flow(graph, set(vin), set(vout), meas_planes) - assert f is not None - assert l_k is not None - pos = vis.get_pos_from_flow(f, l_k) + og = OpenGraph(graph, vin, vout, meas_planes) + causal_flow = og.extract_causal_flow() + pos = vis.get_pos_from_flow(causal_flow) assert pos is not None