From 907b96627ea974830ff9b095e2cc75a1ed12c339 Mon Sep 17 00:00:00 2001 From: Aerylia Date: Thu, 5 Aug 2021 18:23:27 +0300 Subject: [PATCH 1/2] Fixed recursive steiner-gauss bug and added documentation to make the code more legible, plus some refactoring --- pyzx/routing/architecture.py | 28 +++++++++++----- pyzx/routing/steiner.py | 64 ++++++++++++++++++++++++------------ 2 files changed, 62 insertions(+), 30 deletions(-) diff --git a/pyzx/routing/architecture.py b/pyzx/routing/architecture.py index 68294532..e420da19 100644 --- a/pyzx/routing/architecture.py +++ b/pyzx/routing/architecture.py @@ -396,7 +396,8 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): "nodes": nodes, "usable nodes": usable_nodes, "rec nodes": rec_nodes, - "upper": upper + "upper": upper, + "graph_trace": [] } start = self.qubit_map[start] usable_nodes = [self.qubit_map[i] for i in usable_nodes] @@ -407,14 +408,19 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): vertices = [self.vertices[i] for i in usable_nodes] for edge in self.graph.edges(): src, tgt = self.graph.edge_st(edge) + debug_trace["graph_trace"].append((src, tgt)) + # Make sure that the edge is part of the subgraph spanned by vertices if src in vertices and tgt in vertices: - if upper or (src in rec_nodes and tgt in rec_nodes): + # If upper, we don't care about direction + # If one of the vertices in the edge are in rec_nodes, we don't care about direction either. + if upper or (src in rec_nodes or tgt in rec_nodes): distances[(src, tgt)] = (1, [(src, tgt)]) distances[(tgt, src)] = (1, [(tgt, src)]) elif src > tgt: distances[(src, tgt)] = (1, [(src, tgt)]) else: distances[(tgt, src)] = (1, [(tgt, src)]) + # The distance with yourself is 0 for each vertex for v in vertices: distances[(v, v)] = (0, []) for i, v0 in enumerate(vertices+vertices): @@ -426,6 +432,9 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): or distances[(v0, v2)][0] > distances[(v0, v1)][0] + distances[(v1, v2)][0]: distances[(v0, v2)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], distances[(v0, v1)][1] + distances[(v1, v2)][1]) + # If upper, we know that direction doesn't matter, so the reverse is also true. + # this can be found by future iterations so this is a speedup + # Note that when using recursive nodes, we cannot always reverse the direction of the path. if upper: distances[(v2, v0)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], distances[(v2, v1)][1] + distances[(v1, v0)][1]) @@ -437,13 +446,14 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): options = [(node, v, *distances[(v, node)]) for node in nodes for v in (vertices + steiner_pnts) if (v, node) in distances.keys()] if options == []: - print(nodes) - print(steiner_pnts, vertices) - print(distances) - print(usable_nodes) - print(rec_nodes) - print(self.reduce_order) - print(debug_trace) + print("Unable to create a steiner tree. No path can be found.") + print("nodes", nodes) + print("steiner nodes + all vertices", steiner_pnts, vertices) + print("distances", distances) + print("useable nodes", usable_nodes) + print("recursion nodes", rec_nodes) + print("reduce order", self.reduce_order) + print("debug trace", debug_trace) self.visualize("debug.png") print(self.name) best_option = min(options, key=lambda x: x[2]) diff --git a/pyzx/routing/steiner.py b/pyzx/routing/steiner.py index 3097c9d6..b1323427 100644 --- a/pyzx/routing/steiner.py +++ b/pyzx/routing/steiner.py @@ -110,7 +110,8 @@ def steiner_reduce(col, root, nodes, upper): def rec_steiner_gauss(matrix, architecture, full_reduce=False, x=None, y=None, permutation=None, **kwargs): """ - Performs Gaussian elimination that is constraint bij the given architecture + Performs Gaussian elimination that is constraint bij the given architecture according to https://arxiv.org/pdf/1904.00633.pdf + Only works on full rank, square matrices. :param matrix: PyZX Mat2 matrix to be reduced :param architecture: The Architecture object to conform to @@ -136,45 +137,66 @@ def row_add(c0, c1): def steiner_reduce(col, root, nodes, usable_nodes, rec_nodes, upper): generator = steiner_reduce_column(architecture, [row[col] for row in matrix.data], root, nodes, usable_nodes, rec_nodes, upper) cnot = next(generator, None) + tree_nodes = [] while cnot is not None: + tree_nodes.extend(cnot) row_add(*cnot) cnot = next(generator, None) - def rec_step(cols, rows): - size = len(rows) - # Upper triangular form uses the same structure. + return tree_nodes + def rec_step(qubit_removal_order): + # size, p_cols and pivot is needed if the matrix isn't square or of full rank + size = len(qubit_removal_order) p_cols = [] pivot = 0 - rows2 = [r for r in range(len(matrix.data[0])) if r in rows] - cols2 = [c for c in range(len(matrix.data)) if c in cols] + # Order the rows and columns to be ascending. + rows2 = [r for r in range(len(matrix.data[0])) if r in qubit_removal_order] + cols2 = [c for c in range(len(matrix.data)) if c in qubit_removal_order] for i, c in enumerate(cols2): if pivot < size: - nodes = [r for r in rows2[pivot:] if rows2[pivot]==r or matrix.data[r][c] == 1] + # Get all the nodes in the row below the diagonal (rows[i:]) where the entry equals 1 or it is the diagonal + nodes = [r for r in rows2[pivot:] if rows2[pivot]==r or matrix.data[r][c] == 1] + # Perform the steiner_reduce on the current column (c) with root rows2[pivot] and you can use the full matrix steiner_reduce(c, rows2[pivot], nodes, cols2[i:], [], True) if matrix.data[rows2[pivot]][c] == 1: p_cols.append(c) pivot += 1 + #Quick check upper triangular form - should never be printing! + if not all([all([matrix.data[i][j]== 0 for j in range(0, i)]) for i in range(size)]): + print("This should never be printed. If you read this, there is a bug in pyzx/routing/steiner.py line 163") + print() + print("not (upside-down) upper triangular form?") + for row in matrix.data: + print(row) + print("--------") # Full reduce requires the recursion if full_reduce: pivot -= 1 - for i, c in enumerate(cols): + # Get the rows in the same order as the cols, but backwards because we are going backwards + rows = [i for i in reversed(qubit_removal_order) if i in rows2] + for i, c in enumerate(qubit_removal_order): if c in p_cols: - nodes = [r for r in rows if r==rows[pivot] or matrix.data[r][c] == 1] - usable_nodes = cols[i:] - #rec_nodes = list(set([node for edge in architecture.distances["upper"][c][(max(usable_nodes),c)][1] for node in edge])) - #rec_nodes = [n for n in usable_nodes if n in rec_nodes] + # Vertices we can still use steiner trees + usable_nodes = qubit_removal_order[i:] # = R + + # Pick the maximal vertex k in R: k = max(usable_nodes) + # Pick the maximal leaf k' in R: k' = cols[i] = c + # Let W be the set of vertices in the shortest path from k' to k (inclusive) rec_nodes = architecture.shortest_path(c, max(usable_nodes), usable_nodes) - if len(nodes) > 1: - #print("-", c, nodes, cols, rec_nodes) - steiner_reduce(c, rows[pivot], nodes, cols, rec_nodes, False) + # Apply steiner up: + # Pick the nodes of the steiner tree: all rows with a 1 in column c or the pivot row (diagonal) + nodes = [r for r in rows if r==rows[pivot] or matrix.data[r][c] == 1] + + if len(nodes) > 1: # Otherwise it is the diagonal, which is trivial/done + steiner_reduce(c, rows[pivot], nodes, usable_nodes, rec_nodes, False) - # Do recursion on the given nodes. - if len(rec_nodes) > 1: - rec_step(list(reversed(rec_nodes)), rec_nodes) + # Do recursion on the recursive nodes that were allowed to break the the upper triangular form. + if len(rec_nodes) > 1: # Trivial otherwise + rec_step(list(reversed(rec_nodes))) pivot -= 1 - #return rank + # The qubit order is the order in which the spanning tree R is traversed qubit_order = architecture.reduce_order - rec_step(qubit_order, list(reversed(qubit_order))) + rec_step(qubit_order) #print("Done!") def steiner_reduce_column(architecture, col, root, nodes, usable_nodes, rec_nodes, upper): @@ -198,7 +220,7 @@ def steiner_reduce_column(architecture, col, root, nodes, usable_nodes, rec_node else: debug and print("deal with zero root") if next_check is not None and col[next_check[0]] == 0: # root is zero - print("WARNING : Root is 0 => reducing non-pivot column", col) + print("WARNING : Root is 0 => reducing non-pivot column", col, next_check[0]) debug and print("Step 1: remove zeros", col) while next_check is not None: s0, s1 = next_check From 0648d9e0f7f656f335ef5b8cbc7c57a682dd0ccc Mon Sep 17 00:00:00 2001 From: Aerylia Date: Sat, 14 Aug 2021 15:15:06 +0300 Subject: [PATCH 2/2] Added automatic architecture labeling and finding of the reduce order/spanning tree for ease fo use. Included unittests for all architectures. --- pyzx/routing/architecture.py | 483 ++++++++++++++---------------- pyzx/routing/phase_poly.py | 4 +- pyzx/routing/steiner.py | 83 +++-- pyzx/scripts/phase_poly_router.py | 1 - tests/test_cnot_mapper.py | 55 +++- 5 files changed, 304 insertions(+), 322 deletions(-) diff --git a/pyzx/routing/architecture.py b/pyzx/routing/architecture.py index e420da19..fd36a6ae 100644 --- a/pyzx/routing/architecture.py +++ b/pyzx/routing/architecture.py @@ -16,6 +16,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +from math import dist import sys from ..graph.graph import Graph #from pyzx.graph.base import BaseGraph # TODO fix the right graph import - one of many - right backend etc @@ -48,61 +49,109 @@ REC_ARCH, SYCAMORE_LIKE, IBMQ_POUGHKEEPSIE, IBMQ_BOEBLINGEN, IBMQ_SINGAPORE, GOOGLE_SYCAMORE, IBM_ROCHESTER] dynamic_size_architectures = [FULLY_CONNNECTED, LINE, CIRCLE, SQUARE, DENSITY] +hamiltonian_path_architectures = [FULLY_CONNNECTED, LINE, CIRCLE, SQUARE, IBM_QX4, IBM_QX2, IBM_QX3, + IBM_QX5, IBM_Q20_TOKYO, RIGETTI_8Q_AGAVE, RIGETTI_16Q_ASPEN, + IBMQ_POUGHKEEPSIE] -debug = False +debug = True class Architecture(): - def __init__(self, name, coupling_graph=None, coupling_matrix=None, backend=None, qubit_map = None, reduce_order=None, **kwargs): + def __init__(self, name, coupling_graph=None, coupling_matrix=None, backend=None, qubit_map=None, reduce_order=None, **kwargs): """ Class that represents the architecture of the qubits to be taken into account when routing. :param coupling_graph: a PyZX Graph representing the architecture, optional :param coupling_matrix: a 2D numpy array representing the adjacency of the qubits, from which the Graph is created, optional :param backend: The PyZX Graph backend to be used when creating it from the adjacency matrix, optional + :param qubit_map: A qubit placement mapping list such that qubit_map[logical_qubit] = graph_node """ self.name = name if coupling_graph is None: self.graph = Graph(backend=backend) - self._double_edged = True else: self.graph = coupling_graph - self._double_edged = False if coupling_matrix is not None: # build the architecture graph n = coupling_matrix.shape[0] - self.vertices = self.graph.add_vertices(n) - edges = [(self.vertices[row], self.vertices[col]) for row in range(n) for col in range(n) if + self.graph.add_vertices(n) + edges = [(row, col) for row in range(n) for col in range(n) if coupling_matrix[row, col] == 1] self.graph.add_edges(edges) - else: - self.vertices = [v for v in self.graph.vertices()] + + self.vertices = [v for v in self.graph.vertices()] + self.distances = None if qubit_map is not None: self.qubit_map = qubit_map + elif reduce_order is not None: # No qubit map given, but there is a reduce order, so assuming a default i-i mapping. + self.qubit_map = list(range(len(reduce_order))) else: - self.qubit_map = [i for i, v in enumerate(self.vertices)] + self._place_qubits() self.n_qubits = len(self.vertices) self._non_cutting_vertices = None #self.pre_calc_non_cutting_vertices() - self.reduce_order = reduce_order if reduce_order is not None else [i-1 for i in range(self.n_qubits, 0,-1)] + self.reduce_order = self._get_reduce_order() if reduce_order is None else reduce_order self.arities = [(i, len([edge for edge in self.graph.edges() if edge[0] == v])) for i,v in enumerate(self.vertices)] self.arities.sort(key=lambda p: p[1], reverse=True) - self.density = self._density(self.graph.nedges, self.n_qubits) - def pre_calc_distances(self): - self.distances = {"upper": [self.floyd_warshall(until, upper=True) for until, v in enumerate(self.vertices)], - "full": [self.floyd_warshall(until, upper=False) for until, v in enumerate(self.vertices)]} + def qubit2vertex(self, qubit): + # if logical qubit q is stored in vertex v, than self.qubit_map[q] = v + return self.qubit_map[qubit] - def _density(self, edges, qubits): - density = edges/(qubits*(qubits-1)) - if self._double_edged: - return density - return 2*density + def vertex2qubit(self, vertex): + # if logical qubit q is stored in vertex v, than self.qubit_map.index(v) = q + return self.qubit_map.index(vertex) + + def pre_calc_distances(self): + self.distances = {"upper": [self.floyd_warshall(self.vertices[until:], upper=True) for until, v in enumerate(self.vertices)], + "full": [self.floyd_warshall(self.vertices[:until+1], upper=False) for until, v in enumerate(self.vertices)]} + + def _get_reduce_order(self): + #print("qubitmap", self.qubit_map) + vertices = list(sorted(self.vertices, reverse=True, key=self.vertex2qubit)) # sort qubits from large to small + # Pick leaf with largest label every time. + reduce_order = [] + while vertices != []: + all_cutting = self._is_cutting(vertices) # Get which vertices are cutting + # All False indices correspond to a leaf, we want the first one + leaf = all_cutting.index(False) + # Which logical qubit is stored in the leaf node? + reduce_order.append(self.vertex2qubit(vertices[leaf])) + del vertices[leaf] + + return reduce_order + + def _place_qubits(self, start_vertex=None): + # Start at start_vertex and label the graph using DFT and post-ordered labeling + # If no start_vertex is given, start at the node with the highest index + if start_vertex is None: + start_vertex = max(self.vertices) + + visited = [] + # If logical qubit q is stored in physical qubit n, then qubit_map[q]=n + qubit_map = [] + # DFT relabeling + def dft(current): + visited.append(current) # Mark node as visited + # Find the neigbors of current. + neighbors = self.get_neighboring_vertices(current) + #print("current", current) + #print("neighbors", neighbors) + # For all neigbors not yet visited do dft + for n in sorted(neighbors, reverse=True): # Sorting should make it follow the old placements as close as possible + if n not in visited: + dft(n) + # Label current node + qubit_map.append(current) + dft(start_vertex) + #print("visited", visited) + self.qubit_map = qubit_map def pre_calc_non_cutting_vertices(self): # TODO implement this and uncomment line in constructor. - #raise NotImplementedError("pre calculation non cutting vertices") + raise NotImplementedError("pre calculation non cutting vertices") + # TODO refactor using vertex2qubit() and qubit2vertex() and proper naming to make a difference between qubits and vertices qubits = [i for i in range(self.n_qubits)] def collect_non_cutting(qubits): if qubits == []: @@ -117,46 +166,46 @@ def collect_non_cutting(qubits): return all_non_cutting self._non_cutting_vertices = collect_non_cutting(qubits) - def non_cutting_vertices(self, subgraph, pre_calc=False): + def non_cutting_vertices(self, subgraph_vertices, pre_calc=False): # Find the precalculated non-cutting vertices for this subgraph if pre_calc: if self._non_cutting_vertices is None: self.pre_calc_non_cutting_vertices() - if subgraph == []: + if subgraph_vertices == []: return [] subgraphs, non_cutting = zip(*self._non_cutting_vertices) - index = subgraphs.index(subgraph) + index = subgraphs.index(subgraph_vertices) return non_cutting[index] else: if self._non_cutting_vertices is None: self._non_cutting_vertices = {} cur_dict = self._non_cutting_vertices - for q in sorted([v for i, v in enumerate(self.vertices) if i not in subgraph]): + for q in sorted([v for i, v in enumerate(self.vertices) if i not in subgraph_vertices]): if q not in cur_dict.keys(): cur_dict[q] = {} cur_dict = cur_dict[q] if "non_cutting" not in cur_dict.keys(): - cur_dict["non_cutting"] = [subgraph[i] for i, cutting in enumerate(self._is_cutting(qubits=[self.vertices[i] for i in subgraph])) if not cutting] + cur_dict["non_cutting"] = [subgraph_vertices[i] for i, cutting in enumerate(self._is_cutting([self.vertices[i] for i in subgraph_vertices])) if not cutting] #print(self._non_cutting_vertices) return cur_dict["non_cutting"] - def _is_cutting(self, qubits=None): + def _is_cutting(self, vertices=None): # algortihm from https://courses.cs.washington.edu/courses/cse421/04su/slides/artic.pdf and https://www.geeksforgeeks.org/articulation-points-or-cut-vertices-in-a-graph/ # Code written myself. - if qubits is None: - qubits = self.vertices - n_qubits = len(qubits) - discovery_times = [-1]*n_qubits - lows = [None]*n_qubits - index_lookup = {self.vertices[v]:i for i, v in enumerate(qubits)} + if vertices is None: + vertices = self.vertices + number_of_nodes = len(vertices) + discovery_times = [-1]*number_of_nodes + lows = [None]*number_of_nodes + index_lookup = {self.vertices[v]:i for i, v in enumerate(vertices)} self.dfs_counter = 0 - edges = [e for e in self.graph.edges() if e[0] in qubits and e[1] in qubits] + edges = [e for e in self.graph.edges() if e[0] in vertices and e[1] in vertices] edges += [(v2, v1) for v1, v2 in edges] - cutting = [False] * n_qubits - parent = [-1] * n_qubits - def dfs(vertex,): + cutting = [False] * number_of_nodes + parent = [-1] * number_of_nodes + def dfs(vertex): v = index_lookup[vertex] self.dfs_counter += 1 discovery_times[v] = self.dfs_counter @@ -176,7 +225,7 @@ def dfs(vertex,): cutting[v] = True elif v2 != parent[v]: lows[v] = min(lows[v], discovery_times[v2]) - for vertex in qubits: + for vertex in vertices: v = index_lookup[vertex] if discovery_times[v] == -1: dfs(vertex) @@ -187,8 +236,13 @@ def dfs(vertex,): def get_neighboring_qubits(self, qubit): - neighboring_qubits = set([edge[1] for edge in self.graph.edges() if edge[0] == self.vertices[qubit] ] + [edge[0] for edge in self.graph.edges() if edge[1] == self.vertices[qubit]]) - return [self.vertices.index(q) for q in neighboring_qubits] + vertex = self.qubit2vertex(qubit) + return [self.vertex2qubit(q) for q in self.get_neighboring_vertices(vertex)] + + def get_neighboring_vertices(self, vertex): + neighbors = [edge[1] for edge in self.graph.edges() if edge[0] == vertex ] + neighbors += [edge[0] for edge in self.graph.edges() if edge[1] == vertex] + return set(neighbors) def to_quil_device(self): # Only required here @@ -211,51 +265,53 @@ def visualize(self, filename=None): filename = self.name + ".png" plt.savefig(filename) - def floyd_warshall(self, exclude_excl, upper=True): + def floyd_warshall(self, subgraph_vertices=None, upper=True, rec_vertices=[]): """ Implementation of the Floyd-Warshall algorithm to calculate the all-pair distances in a given graph - :param exclude_excl: index up to which qubit should be excluded from the distances + :param subgraph_vertices: subset of vertices to calculate the distances over :param upper: whether use bidirectional edges or only ordered edges (src, tgt) such that src > tgt, default True + :param rec_vertices: A list of vertices for which you are allowed out-of-order edges if not upper :return: a dict with for each pair of qubits in the graph, a tuple with their distance and the corresponding shortest path """ # https://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm distances = {} - vertices = self.vertices[exclude_excl:] if upper else self.vertices[:exclude_excl + 1] + vertices = subgraph_vertices if subgraph_vertices is not None else self.vertices for edge in self.graph.edges(): src, tgt = self.graph.edge_st(edge) if src in vertices and tgt in vertices: - if upper: + if upper or (src in rec_vertices and tgt in rec_vertices): distances[(src, tgt)] = (1, [(src, tgt)]) distances[(tgt, src)] = (1, [(tgt, src)]) - elif src > tgt: + elif self.vertex2qubit(src) > self.vertex2qubit(tgt): distances[(src, tgt)] = (1, [(src, tgt)]) else: distances[(tgt, src)] = (1, [(tgt, src)]) for v in vertices: distances[(v, v)] = (0, []) - for i, v0 in enumerate(vertices): - for j, v1 in enumerate(vertices if upper else vertices[:i + 1]): - for v2 in vertices if upper else vertices[: i + j + 1]: - if (v0, v1) in distances.keys(): - if (v1, v2) in distances.keys(): - if (v0, v2) not in distances.keys() or distances[(v0, v2)][0] > distances[(v0, v1)][0] + \ - distances[(v1, v2)][0]: - distances[(v0, v2)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], - distances[(v0, v1)][1] + distances[(v1, v2)][1]) - if upper: - distances[(v2, v0)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], - distances[(v2, v1)][1] + distances[(v1, v0)][1]) + for v0 in vertices+vertices: + for v1 in vertices: + for v2 in vertices: + connecting_path_found = (v0, v1) in distances.keys() and (v1, v2) in distances.keys() + new_path_found = connecting_path_found and (v0, v2) not in distances.keys() + shorter_path_found = connecting_path_found and not new_path_found and distances[(v0, v2)][0] > distances[(v0, v1)][0] + distances[(v1, v2)][0] + # If I can go from v0 to v1 and from v1 to v2, I can go from v0 to v2. => I can always take this path. + if new_path_found or shorter_path_found: + distances[(v0, v2)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], + distances[(v0, v1)][1] + distances[(v1, v2)][1]) + if upper: + distances[(v2, v0)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], + distances[(v2, v1)][1] + distances[(v1, v0)][1]) return distances - def shortest_path(self, start, end, subgraph=None): + def shortest_path(self, start_qubit, end_qubit, qubits_to_use=None): #print("shortest path") - if subgraph is None: + if qubits_to_use is None: nodes = self.vertices else: - nodes = [self.vertices[n] for n in subgraph] - start = self.vertices[start] - end = self.vertices[end] + nodes = [self.qubit2vertex(n) for n in qubits_to_use] + start = self.qubit2vertex(start_qubit) + end = self.qubit2vertex(end_qubit) queue = [(start, [start])] visited = [start] @@ -265,33 +321,36 @@ def shortest_path(self, start, end, subgraph=None): node, path = queue.pop() if node == end: - #print(start, end, subgraph) + #print(start, end, qubits_to_use) #self.visualize("debug.png") #print(path) return path edges = [edge for edge in self.graph.edges() if node in edge] neighbors = [n for edge in edges for n in edge if n != node and n in nodes] #print(node, neighbors) - for node in neighbors: - if node not in visited: - queue.append((node, path + [node])) - visited.append(node) + for new_node in neighbors: + if new_node not in visited: + queue.append((new_node, path + [new_node])) + visited.append(new_node) debug_trace.append(queue[-1]) - print(start, end, subgraph) + print("No shortest path found!") + print(start_qubit, end_qubit, qubits_to_use) + print(start, end, nodes) print(*self.graph.edges()) + print([self.qubit2vertex(v) for v in self.reduce_order]) print("Trace:", *debug_trace, sep="\n") self.visualize("debug.png") print(self.name) return None - def steiner_tree(self, start, nodes, upper=True): + def steiner_tree(self, start_qubit, qubits_to_use, upper=True): """ Approximates the steiner tree given the architecture, a root qubit and the other qubits that should be present. This is done using the pre-calculated all-pairs shortest distance and Prim's algorithm for creating a minimum spanning tree - :param start: The index of the root qubit to be used - :param nodes: The indices of the other qubits that should be present in the steiner tree + :param start_qubit: The index of the root qubit to be used + :param qubits_to_use: The indices of the other qubits that should be present in the steiner tree :param upper: Whether the steiner tree is used for creating an upper triangular matrix or a full reduction. - :yields: First yields all edges from the tree top-to-bottom, finished with None, then yields all edges from the tree bottom-up, finished with None. + :yields: First yields all edges from the tree top-to-bottom, finished with None, then yields all edges from the tree bottom-up, finished with None. Where the edges are qubit indices, not vertex indices """ # Approximated by calculating the all-pairs shortest paths and then solving the mininum spanning tree over the subset of vertices and their respective shortest paths. # https://en.wikipedia.org/wiki/Steiner_tree_problem#Approximating_the_Steiner_tree @@ -302,9 +361,10 @@ def steiner_tree(self, start, nodes, upper=True): # returns an iterator that walks the steiner tree, yielding (adj_node, leaf) pairs. If the walk is finished, it yields None if self.distances is None: self.pre_calc_distances() + start = self.qubit2vertex(start_qubit) + nodes = [self.qubit2vertex(q) for q in qubits_to_use] state = [start, [n for n in nodes]] root = start - # TODO deal with qubit mapping vertices = [root] edges = [] debug and print(root, upper, nodes) @@ -331,38 +391,44 @@ def steiner_tree(self, start, nodes, upper=True): vs = {root} n_edges = len(edges) yielded_edges = set() - debug_count = 0 - yield_count = 0 - warning = 0 + if debug: + debug_count = 0 + yield_count = 0 + warning = 0 + while len(yielded_edges) < n_edges: es = [e for e in edges for v in vs if e[0] == v] old_vs = [v for v in vs] yielded = False for edge in es: - yield edge + yield (self.vertex2qubit(edge[0]), self.vertex2qubit(edge[1])) vs.add(edge[1]) if edge in yielded_edges: print("DOUBLE yielding! - should not be possible!") yielded_edges.add(edge) - yielded = True - yield_count += 1 + if debug: + yielded = True + yield_count += 1 + [vs.remove(v) for v in old_vs] - if not yielded: - debug and print("leaf!") - debug_count += 1 - if debug_count > len(vertices): - print("infinite loop!", warning) + if debug: + if not yielded: + debug and print("leaf!") + debug_count += 1 + if debug_count > len(vertices): + print("infinite loop!", warning) + warning += 1 + if yield_count > len(edges): + print("Yielded more edges than existing... This should not be possible!", warning) warning += 1 - if yield_count > len(edges): - print("Yielded more edges than existing... This should not be possible!", warning) - warning += 1 - if warning > 5: - print(state, yielded_edges) - # input("note it down") - break + if warning > 5: + print(state, yielded_edges) + # input("note it down") + break yield None # Walk the tree bottom up to remove all ones. - yield_count = 0 + if debug: + yield_count = 0 while len(edges) > 0: # find leaf nodes: debug and print(vertices, steiner_pnts, edges) @@ -372,72 +438,47 @@ def steiner_tree(self, start, nodes, upper=True): for v in vs_to_consider: # Get the edge that is connected to this leaf node for edge in [e for e in edges if e[1] == v]: - yield edge + yield (self.vertex2qubit(edge[0]), self.vertex2qubit(edge[1])) edges.remove(edge) - yielded = True - yield_count += 1 - # yield map(lambda i: self.qubit_map[i], edge) - if not yielded: - print("Infinite loop!", warning) - warning += 1 - if yield_count > n_edges: - print("Yielded more edges than existing again... This should not be possible!!", warning) - warning += 1 - if warning > 10: - print(state, edges, yield_count) - # input("Note it down!") - break + if debug: + yielded = True + yield_count += 1 + if debug: + if not yielded: + print("Infinite loop!", warning) + warning += 1 + if yield_count > n_edges: + print("Yielded more edges than existing again... This should not be possible!!", warning) + warning += 1 + if warning > 10: + print(state, edges, yield_count) + # input("Note it down!") + break yield None - def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): + def rec_steiner_tree(self, start_qubit, terminal_qubits, usable_qubits, rec_qubits, upper=True): + if not all([q in usable_qubits for q in terminal_qubits]): + raise Exception("Terminals not in the subgraph") # Builds the steiner tree with start as root, contains at least nodes and at most useable_nodes debug_trace = { - "start":start, - "nodes": nodes, - "usable nodes": usable_nodes, - "rec nodes": rec_nodes, + "start":start_qubit, + "terminal qubits": terminal_qubits, + "usable qubits": usable_qubits, + "rec qubits": rec_qubits, "upper": upper, "graph_trace": [] } - start = self.qubit_map[start] - usable_nodes = [self.qubit_map[i] for i in usable_nodes] - nodes = [self.qubit_map[i] for i in nodes] - rec_nodes = [self.qubit_map[i] for i in rec_nodes] + start = self.qubit2vertex(start_qubit) + usable_nodes = [self.qubit2vertex(i) for i in usable_qubits] + nodes = [self.qubit2vertex(i) for i in terminal_qubits] + rec_nodes = [self.qubit2vertex(i) for i in rec_qubits] # Calculate all-pairs shortest path - distances = {} - vertices = [self.vertices[i] for i in usable_nodes] - for edge in self.graph.edges(): - src, tgt = self.graph.edge_st(edge) - debug_trace["graph_trace"].append((src, tgt)) - # Make sure that the edge is part of the subgraph spanned by vertices - if src in vertices and tgt in vertices: - # If upper, we don't care about direction - # If one of the vertices in the edge are in rec_nodes, we don't care about direction either. - if upper or (src in rec_nodes or tgt in rec_nodes): - distances[(src, tgt)] = (1, [(src, tgt)]) - distances[(tgt, src)] = (1, [(tgt, src)]) - elif src > tgt: - distances[(src, tgt)] = (1, [(src, tgt)]) - else: - distances[(tgt, src)] = (1, [(tgt, src)]) - # The distance with yourself is 0 for each vertex - for v in vertices: - distances[(v, v)] = (0, []) - for i, v0 in enumerate(vertices+vertices): - for j, v1 in enumerate(vertices): - for v2 in vertices: - if (v0, v1) in distances.keys(): - if (v1, v2) in distances.keys(): - if (v0, v2) not in distances.keys() \ - or distances[(v0, v2)][0] > distances[(v0, v1)][0] + distances[(v1, v2)][0]: - distances[(v0, v2)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], - distances[(v0, v1)][1] + distances[(v1, v2)][1]) - # If upper, we know that direction doesn't matter, so the reverse is also true. - # this can be found by future iterations so this is a speedup - # Note that when using recursive nodes, we cannot always reverse the direction of the path. - if upper: - distances[(v2, v0)] = (distances[(v0, v1)][0] + distances[(v1, v2)][0], - distances[(v2, v1)][1] + distances[(v1, v0)][1]) + distances = self.floyd_warshall(usable_nodes, upper=upper, rec_vertices=rec_nodes) + if debug: + debug_trace["root vtx"] = start + debug_trace["terminals"] = nodes + debug_trace["subgraph"] = usable_nodes + debug_trace["recursion nodes"] = rec_nodes # Build the spanning tree of shortest paths with root start, containing at least nodes vertices = [start] edges = [] @@ -447,15 +488,21 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): (v, node) in distances.keys()] if options == []: print("Unable to create a steiner tree. No path can be found.") - print("nodes", nodes) - print("steiner nodes + all vertices", steiner_pnts, vertices) - print("distances", distances) - print("useable nodes", usable_nodes) - print("recursion nodes", rec_nodes) - print("reduce order", self.reduce_order) - print("debug trace", debug_trace) - self.visualize("debug.png") - print(self.name) + if debug: + print("vertices", nodes) + print("steiner nodes + all vertices", steiner_pnts, vertices) + print("distances:") + for pair, distance in distances.items(): + print("\t", pair, distance[0], "\t", distance[1]) + + print("Edges", [e for e in self.graph.edges()]) + print("useable vertices", usable_nodes) + print("recursion vertices", rec_nodes) + print("reduce order", self.reduce_order) + print("qubit map", self.qubit_map) + print("debug trace", debug_trace) + self.visualize("debug.png") + print(self.name) best_option = min(options, key=lambda x: x[2]) vertices.append(best_option[0]) edges += best_option[3] @@ -471,7 +518,7 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): es = [e for e in edges for v in vs if e[0] == v] # Find all vertices connected to previously yielded vertices old_vs = [v for v in vs] for edge in es: # yield the corresponding edges. - yield (self.qubit_map.index(edge[0]), self.qubit_map.index(edge[1])) + yield (self.vertex2qubit(edge[0]), self.vertex2qubit(edge[1])) vs.add(edge[1]) yielded_edges.add(edge) [vs.remove(v) for v in old_vs] @@ -484,7 +531,7 @@ def rec_steiner_tree(self, start, nodes, usable_nodes, rec_nodes, upper=True): for v in vs_to_consider: # Get the edge that is connected to this leaf node for edge in [e for e in edges if e[1] == v]: - yield (self.qubit_map.index(edge[0]), self.qubit_map.index(edge[1])) # yield it + yield (self.vertex2qubit(edge[0]), self.vertex2qubit(edge[1])) # yield it edges.remove(edge) # Remove it from the steiner tree yield None # Signal done @@ -669,8 +716,8 @@ def create_ibmq_singapore(backend=None, name=None, **kwargs): graph.add_edges(edges) reduce_order = [19, 18, 17, 15, 16, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 3, 4, 2, 1, 0] if name is not None: - return Architecture(name=IBMQ_BOEBLINGEN, coupling_graph=graph, backend=backend, reduce_order=reduce_order, **kwargs) - return Architecture(name=IBMQ_SINGAPORE, coupling_graph=graph, backend=backend, reduce_order=reduce_order, **kwargs) + return Architecture(name=IBMQ_BOEBLINGEN, coupling_graph=graph, backend=backend, **kwargs) + return Architecture(name=IBMQ_SINGAPORE, coupling_graph=graph, backend=backend, **kwargs) def create_rigetti_19q_acorn_architecture(backend=None, **kwargs): graph = Graph(backend=backend) @@ -679,7 +726,7 @@ def create_rigetti_19q_acorn_architecture(backend=None, **kwargs): extra_edges = [(8,9), (0,18), (2,16), (4,14), (6, 12)] edges += [(vertices[v1], vertices[v2]) for v1, v2 in extra_edges] graph.add_edges(edges) - return Architecture(RIGETTI_19Q_ACORN, coupling_graph=graph, reduce_order=[19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 8, 9, 7, 6, 5, 4, 3, 2, 1, 0], backend=backend, **kwargs) + return Architecture(RIGETTI_19Q_ACORN, coupling_graph=graph, backend=backend, **kwargs) def create_rigetti_16q_aspen_architecture(backend=None, **kwargs): @@ -699,7 +746,7 @@ def create_sycamore_like(backend=None, **kwargs): extra_edges = [(1,2),(2,6), (4,5), (4,14),(5,12), (6,11), (7,10), (13,14), (12,16), (11,18), (10,17), (17,18)] edges += [(vertices[v1], vertices[v2]) for v1, v2 in extra_edges] graph.add_edges(edges) - return Architecture(SYCAMORE_LIKE, coupling_graph=graph, backend=backend, reduce_order=[19, 17, 18, 16,15,13,14,12,11,10,9,8,7,6,4,5,3,1,2,0], **kwargs) + return Architecture(SYCAMORE_LIKE, coupling_graph=graph, backend=backend, **kwargs) def create_rigetti_8q_agave_architecture(**kwargs): m = np.array([ @@ -727,7 +774,8 @@ def create_recursive_architecture(**kwargs): [0, 0, 1, 0, 0, 1, 1, 0, 1],#7 [0, 0, 0, 0, 0, 0, 0, 1, 0] #8 ]) - return Architecture(name=REC_ARCH, coupling_matrix=m, reduce_order=[8,6,4,3,5,7,1,2,0], **kwargs) + #reduce_order=[8,6,4,3,5,7,1,2,0] + return Architecture(name=REC_ARCH, coupling_matrix=m, **kwargs) def create_fully_connected_architecture(n_qubits=None, **kwargs): if n_qubits is None: @@ -754,8 +802,6 @@ def create_dynamic_density_hamiltonian_architecture(n_qubits, density_prob=0.1, def create_dynamic_density_tree_architecture(n_qubits, density_prob=0.1, backend=None, **kwargs): import networkx as nx - #graph = Graph(backend=backend) - #vertices = graph.add_vertices(n_qubits) graph = nx.Graph() vertices = [i for i in range(n_qubits)] graph.add_nodes_from(vertices) @@ -790,32 +836,10 @@ def create_dynamic_density_tree_architecture(n_qubits, density_prob=0.1, backend edges.extend([possible_edges[i] for i in indices]) #print(*sorted([(e1,e2) if e1 > e2 else (e2,e1) for e1,e2 in edges], key=lambda p: p[0]), sep="\n") graph.add_edges_from(edges) - - # Number the edges - numbered = {} - numbered = {v:i for i,v in enumerate(nx.dfs_postorder_nodes(graph, source=root))} # Make the coupling graph and adjust the numbering m = nx.to_numpy_array(graph) - perm = [numbered[i] for i in range(n_qubits)] - perm = [perm.index(i) for i in range(n_qubits)] - m = m[perm][:, perm] - spanning_tree = nx.dfs_tree(graph, source=root) - #print(*[[numbered[n] for n in edge] for edge in spanning_tree.edges()]) - g = Graph(backend=backend) - vertices = g.add_vertices(n_qubits) - g.add_edges([(vertices[v1], vertices[v2]) for v1, v2 in [[numbered[n] for n in edge] for edge in spanning_tree.edges()]]) - arch = Architecture(name="temp", coupling_graph=g, backend=backend, **kwargs) - subgraph = [i for i in range(n_qubits)] - reduce_order = [] - while subgraph != []: - non_cutting = arch.non_cutting_vertices(subgraph) - node = max(non_cutting) - reduce_order.append(node) - subgraph.remove(node) - #arch.reduce_order = reduce_order - #reduce_order = [] name = dynamic_size_architecture_name(DENSITY+str(density_prob), n_qubits) - arch = Architecture(name, coupling_matrix=m, reduce_order=reduce_order, **kwargs) + arch = Architecture(name, coupling_matrix=m, **kwargs) return arch def create_dynamic_density_architecture(n_qubits, density_prob=0.1, backend=None, **kwargs): @@ -831,58 +855,9 @@ def create_dynamic_density_architecture(n_qubits, density_prob=0.1, backend=None graph.add_edge(*disconnected[to_add]) connectivity = nx.all_pairs_node_connectivity(graph) disconnected = [(v1, v2) for v1,d in connectivity.items() for v2, score in d.items() if score == 0] - # Number the edges - numbered = {} - numbered = {v:i for i,v in enumerate(nx.dfs_postorder_nodes(graph))} - # Make the coupling graph and adjust the numbering m = nx.to_numpy_array(graph) - perm = [numbered[i] for i in range(n_qubits)] - perm = [perm.index(i) for i in range(n_qubits)] - m = m[perm][:, perm] - # Find the reduce order - """ - g = nx.convert_matrix.from_numpy_array(m) - reduce_order = [] - successors = nx.dfs_successors(g, source=0) - def ordered_dfs(node): - if node not in successors.keys(): # leaf node - return [node] - else: - children = successors[node] - return sum([ordered_dfs(child) for child in sorted(children, reverse=True)], []) + [node] - reduce_order = ordered_dfs(0) - """ - spanning_tree = nx.dfs_tree(graph) - #print(*[[numbered[n] for n in edge] for edge in spanning_tree.edges()]) - g = Graph(backend=backend) - vertices = g.add_vertices(n_qubits) - g.add_edges([(vertices[v1], vertices[v2]) for v1, v2 in [[numbered[n] for n in edge] for edge in spanning_tree.edges()]]) - arch = Architecture(name="temp", coupling_graph=g, backend=backend, **kwargs) - subgraph = [i for i in range(n_qubits)] - reduce_order = [] - while subgraph != []: - non_cutting = arch.non_cutting_vertices(subgraph) - node = max(non_cutting) - reduce_order.append(node) - subgraph.remove(node) - #arch.reduce_order = reduce_order - #reduce_order = [] name = dynamic_size_architecture_name(DENSITY+str(density_prob), n_qubits) - arch = Architecture(name, coupling_matrix=m, reduce_order=reduce_order, **kwargs) - """ - subgraph = [i for i in range(n_qubits)] - reduce_order = [] - while subgraph != []: - non_cutting = arch.non_cutting_vertices(subgraph) - node = max(non_cutting) - reduce_order.append(node) - subgraph.remove(node) - arch.reduce_order = reduce_order - """ - #arch.visualize("dynamic.png") - #print(reduce_order) - #exit(42) - #return Architecture(name, coupling_matrix=m, reduce_order=reduce_order, **kwargs) + arch = Architecture(name, coupling_matrix=m, **kwargs) return arch def create_ibm_rochester(backend=None, **kwargs): @@ -892,15 +867,7 @@ def create_ibm_rochester(backend=None, **kwargs): extra_edges = [(7, 8), (7,20), (14,15), (14,44), (17,18), (17, 28), (0,22), (30,31), (30,42), (37, 38), (40,41), (40,51)] edges += [(vertices[v1], vertices[v2]) for v1, v2 in extra_edges] graph.add_edges(edges) - reduce_order = [] - i = 52 - for j in reversed([7, 14, 17,30,37, 40]): - reduce_order += list(range(i, j+1, -1)) - reduce_order += [j, j+1] - i = j-1 - reduce_order += list(range(i,-1, -1)) - print(reduce_order) - return Architecture(IBM_ROCHESTER, coupling_graph=graph, reduce_order=reduce_order, backend=backend, **kwargs) + return Architecture(IBM_ROCHESTER, coupling_graph=graph, backend=backend, **kwargs) def create_google_sycamore(backend=None, **kwargs): graph = Graph(backend=backend) @@ -917,23 +884,7 @@ def create_google_sycamore(backend=None, **kwargs): extra_edges += list(zip([1,6,21,33,47,50,12,32,31],[4,7,24,38, 48,51,13,43,33])) edges += [(vertices[v1], vertices[v2]) for v1, v2 in extra_edges] graph.add_edges(edges) - reduce_order = [] - i = 52 - for j in reversed([6,12,32, 47,50]): - if j == 32: - reduce_order += list(range(i, j+1, -1)) - reduce_order += [j,31, j+1] - i = j-2 - else: - reduce_order += list(range(i, j+1, -1)) - reduce_order += [j, j+1] - i = j-1 - reduce_order += list(range(i,-1, -1)) - #print(reduce_order) - #print(len(edges)) - #print(*edges, sep='\n') - arch = Architecture(GOOGLE_SYCAMORE, coupling_graph=graph, reduce_order=reduce_order, backend=backend, **kwargs) - #arch.visualize() + arch = Architecture(GOOGLE_SYCAMORE, coupling_graph=graph, backend=backend, **kwargs) return arch def create_architecture(name, **kwargs): @@ -947,7 +898,7 @@ def create_architecture(name, **kwargs): arch_dict[IBMQ_POUGHKEEPSIE] = create_ibmq_poughkeepsie arch_dict[IBMQ_SINGAPORE] = create_ibmq_singapore arch_dict[IBMQ_BOEBLINGEN] = lambda **kwargs : create_ibmq_singapore(name=IBMQ_BOEBLINGEN, **kwargs) - arch_dict[DENSITY] = create_dynamic_density_tree_architecture + arch_dict[DENSITY] = create_dynamic_density_tree_architecture # TODO maybe use a different one arch_dict[GOOGLE_SYCAMORE] = create_google_sycamore arch_dict[IBM_ROCHESTER] = create_ibm_rochester if name == SQUARE: diff --git a/pyzx/routing/phase_poly.py b/pyzx/routing/phase_poly.py index ffa84371..36e3bc16 100644 --- a/pyzx/routing/phase_poly.py +++ b/pyzx/routing/phase_poly.py @@ -692,7 +692,9 @@ def place_cnot(control, target): def base_recurse(cols_to_use, qubits_to_use): if qubits_to_use != [] and cols_to_use != []: # Select edge qubits - qubits = architecture.non_cutting_vertices(qubits_to_use) + vertices_to_use = [architecture.qubit2vertex(q) for q in qubits_to_use] + vertices = architecture.non_cutting_vertices(vertices_to_use) + qubits = [architecture.vertex2qubit(v) for v in vertices] # Pick the qubit where the recursion split will be most skewed. if zeroes_rec: diff --git a/pyzx/routing/steiner.py b/pyzx/routing/steiner.py index b1323427..2d62be91 100644 --- a/pyzx/routing/steiner.py +++ b/pyzx/routing/steiner.py @@ -124,80 +124,75 @@ def rec_steiner_gauss(matrix, architecture, full_reduce=False, x=None, y=None, p permutation = [i for i in range(len(matrix.data))] else: matrix = Mat2([[row[i] for i in permutation] for row in matrix.data]) - #print(matrix) + def row_add(c0, c1): matrix.row_add(c0, c1) #c0 = permutation[c0] #c1 = permutation[c1] debug and print("Reducing", c0, c1) - c0 = architecture.qubit_map[c0] - c1 = architecture.qubit_map[c1] + #c0 = architecture.qubit2vertex(c0) + #c1 = architecture.qubit2vertex(c1) if x != None: x.row_add(c0, c1) if y != None: y.col_add(c1, c0) + def steiner_reduce(col, root, nodes, usable_nodes, rec_nodes, upper): + if not all([q in usable_nodes for q in nodes]): + raise Exception("Terminals not in the subgraph "+ str(upper) + str((col, root, nodes, usable_nodes, rec_nodes))+ "\n"+str(debug_trace)) generator = steiner_reduce_column(architecture, [row[col] for row in matrix.data], root, nodes, usable_nodes, rec_nodes, upper) cnot = next(generator, None) tree_nodes = [] while cnot is not None: + if cnot[0] not in usable_nodes+rec_nodes or cnot[1] not in usable_nodes + rec_nodes: + raise Exception("Steiner tree not sticking to constraints") tree_nodes.extend(cnot) row_add(*cnot) cnot = next(generator, None) return tree_nodes + def rec_step(qubit_removal_order): # size, p_cols and pivot is needed if the matrix isn't square or of full rank size = len(qubit_removal_order) - p_cols = [] - pivot = 0 # Order the rows and columns to be ascending. - rows2 = [r for r in range(len(matrix.data[0])) if r in qubit_removal_order] - cols2 = [c for c in range(len(matrix.data)) if c in qubit_removal_order] - for i, c in enumerate(cols2): - if pivot < size: - # Get all the nodes in the row below the diagonal (rows[i:]) where the entry equals 1 or it is the diagonal - nodes = [r for r in rows2[pivot:] if rows2[pivot]==r or matrix.data[r][c] == 1] - # Perform the steiner_reduce on the current column (c) with root rows2[pivot] and you can use the full matrix - steiner_reduce(c, rows2[pivot], nodes, cols2[i:], [], True) - if matrix.data[rows2[pivot]][c] == 1: - p_cols.append(c) - pivot += 1 + pivots = list(sorted(qubit_removal_order)) + for pivot_idx, c in enumerate(pivots): + # Get all the nodes in the row below the diagonal (rows[i:]) where the entry equals 1 or it is the diagonal + nodes = [r for r in pivots[pivot_idx:] if c==r or matrix.data[r][c] == 1] + # Perform the steiner_reduce on the current column (c) with root rows2[pivot] and you can use the full matrix + steiner_reduce(c, c, nodes, pivots[pivot_idx:], [], True) #Quick check upper triangular form - should never be printing! if not all([all([matrix.data[i][j]== 0 for j in range(0, i)]) for i in range(size)]): - print("This should never be printed. If you read this, there is a bug in pyzx/routing/steiner.py line 163") + print("This should never be printed. If you read this, there is a bug around pyzx/routing/steiner.py line 163") print() - print("not (upside-down) upper triangular form?") + print("not upper triangular form?") for row in matrix.data: print(row) print("--------") # Full reduce requires the recursion if full_reduce: - pivot -= 1 - # Get the rows in the same order as the cols, but backwards because we are going backwards - rows = [i for i in reversed(qubit_removal_order) if i in rows2] + # We precalculated the maximal leafs in R (in the qubit_removal_order) for i, c in enumerate(qubit_removal_order): - if c in p_cols: - # Vertices we can still use steiner trees - usable_nodes = qubit_removal_order[i:] # = R - - # Pick the maximal vertex k in R: k = max(usable_nodes) - # Pick the maximal leaf k' in R: k' = cols[i] = c - # Let W be the set of vertices in the shortest path from k' to k (inclusive) - rec_nodes = architecture.shortest_path(c, max(usable_nodes), usable_nodes) - - # Apply steiner up: - # Pick the nodes of the steiner tree: all rows with a 1 in column c or the pivot row (diagonal) - nodes = [r for r in rows if r==rows[pivot] or matrix.data[r][c] == 1] - - if len(nodes) > 1: # Otherwise it is the diagonal, which is trivial/done - steiner_reduce(c, rows[pivot], nodes, usable_nodes, rec_nodes, False) - - # Do recursion on the recursive nodes that were allowed to break the the upper triangular form. - if len(rec_nodes) > 1: # Trivial otherwise - rec_step(list(reversed(rec_nodes))) - pivot -= 1 + # Vertices we can still use steiner trees + usable_nodes = qubit_removal_order[i:] # = R + + # Pick the maximal vertex k in R: k = max(usable_nodes) + # Pick the maximal leaf k' in R: k' = cols[i] = c + # Let W be the set of vertices in the shortest path from k' to k (inclusive) + path = architecture.shortest_path(c, max(usable_nodes))#, usable_nodes) + rec_qubits = [architecture.vertex2qubit(v) for v in path] + + # Apply steiner up: + # Pick the nodes of the steiner tree: all rows with a 1 in column c or the pivot row (diagonal) + nodes = [r for r in usable_nodes if (r==c or matrix.data[r][c] == 1)] + + if len(nodes) > 1: # Otherwise it is the diagonal, which is trivial/done + steiner_reduce(c, c, nodes, usable_nodes, rec_qubits, False) + + # Do recursion on the recursive nodes that were allowed to break the the upper triangular form. + if len(rec_qubits) > 1: # Trivial otherwise + rec_step(list(reversed(rec_qubits)))#[n for n in rows if n in rec_qubits]) + # The qubit order is the order in which the spanning tree R is traversed - qubit_order = architecture.reduce_order - rec_step(qubit_order) - #print("Done!") + rec_step(architecture.reduce_order) def steiner_reduce_column(architecture, col, root, nodes, usable_nodes, rec_nodes, upper): steiner_tree = architecture.rec_steiner_tree(root, nodes, usable_nodes, rec_nodes, upper) diff --git a/pyzx/scripts/phase_poly_router.py b/pyzx/scripts/phase_poly_router.py index 52e1b452..c5413c83 100644 --- a/pyzx/scripts/phase_poly_router.py +++ b/pyzx/scripts/phase_poly_router.py @@ -192,7 +192,6 @@ def main(args): results_df["matroid"] = method results_df["root_heuristic"] = root_heuristic results_df["split_heuristic"] = split_heuristic - results_df["arch_density"] = architecture.density results_df["architecture"] = architecture.name results_df.set_index(["idx", "mode", "file"], inplace=True, append=True) all_results.append(results_df) diff --git a/tests/test_cnot_mapper.py b/tests/test_cnot_mapper.py index dee91083..152f449d 100644 --- a/tests/test_cnot_mapper.py +++ b/tests/test_cnot_mapper.py @@ -8,7 +8,7 @@ from pyzx.linalg import Mat2 from pyzx.routing.cnot_mapper import gauss, STEINER_MODE, GAUSS_MODE, GENETIC_STEINER_MODE, PSO_STEINER_MODE, cnot_fitness_func, sequential_gauss -from pyzx.routing.architecture import create_architecture, REC_ARCH, SQUARE, IBMQ_SINGAPORE +from pyzx.routing.architecture import Architecture, create_architecture, REC_ARCH, SQUARE, IBMQ_SINGAPORE, RIGETTI_16Q_ASPEN, architectures, dynamic_size_architectures, DENSITY, hamiltonian_path_architectures from pyzx.parity_maps import CNOT_tracker, build_random_parity_map from pyzx.machine_learning import GeneticAlgorithm from pyzx.circuit import CNOT @@ -19,8 +19,8 @@ class TestSteiner(unittest.TestCase): def setUp(self): - self.n_tests = 5 - self.arch = create_architecture(IBMQ_SINGAPORE) + self.n_tests = 10 + self.arch = create_architecture(SQUARE, n_qubits=9) #Needs to have a square number of qubits to test the square architecture. self.n_qubits = self.arch.n_qubits depth = 20 self.circuit = [CNOT_tracker(self.arch.n_qubits) for _ in range(self.n_tests)] @@ -31,11 +31,17 @@ def setUp(self): for g in c.gates: self.aggr_circ.add_gate(g) - def assertGates(self, circuit): + def assertGates(self, circuit, architecture=None): + if architecture is None: + architecture = self.arch for gate in circuit.gates: if hasattr(gate, "name") and gate.name == "CNOT": - edge = (gate.control, gate.target) - self.assertTrue(edge in self.arch.graph.edges() or tuple(reversed(edge)) in self.arch.graph.edges()) + qubits = (gate.control, gate.target) + edge = (architecture.qubit2vertex(qubits[0]), architecture.qubit2vertex(qubits[1])) + edges = list(architecture.graph.edges()) + edges += [tuple(reversed(edge)) for edge in edges] + self.assertIn(edge, edges) + #self.assertTrue(edge in architecture.graph.edges() or tuple(reversed(edge)) in architecture.graph.edges()) def assertMat2Equal(self, m1, m2, triangle=False): if triangle: @@ -89,6 +95,33 @@ def test_upper_shortest_path(self): self.assertEqual(distance, len(path)) """ + def test_architectures(self): + for i in range(self.n_tests): + for name in architectures: + #print("Testing", name) + if name in dynamic_size_architectures: + if name not in hamiltonian_path_architectures: + #continue + # This is DENSITY and it gets in an infinite loop when building steiner trees! + architecture = create_architecture(name, n_qubits=self.n_qubits) + else: + try: + architecture = create_architecture(name, n_qubits=self.n_qubits) + except KeyError as e: + if name == SQUARE: + print("WARNING! skipping SQUARE architecture due to non-square number of qubits:", self.n_qubits) + else: + raise e + elif name in hamiltonian_path_architectures: + architecture = create_architecture(name) + else: + architecture = create_architecture(name) + architecture.visualize(name+".png") + if self.n_qubits == architecture.n_qubits: + with self.subTest(i=i, arch=architecture.name): + self.do_gauss(STEINER_MODE, self.matrix[i], architecture=architecture) + + def test_all_cnots_valid(self): for i in range(self.n_tests): with self.subTest(i=i): @@ -97,11 +130,13 @@ def test_all_cnots_valid(self): _ = gauss(STEINER_MODE, matrix, architecture=self.arch, full_reduce=True, y=circuit) self.assertGates(circuit) - def do_gauss(self, mode, array, full_reduce=True, with_assert=True): + def do_gauss(self, mode, array, full_reduce=True, with_assert=True, architecture=None): + if architecture is None: + architecture = self.arch circuit = CNOT_tracker(self.arch.n_qubits) matrix = Mat2(np.copy(array)) - rank = gauss(mode, matrix, architecture=self.arch, full_reduce=full_reduce, y=circuit) - with_assert and mode == STEINER_MODE and self.assertGates(circuit) + rank = gauss(mode, matrix, architecture=architecture, full_reduce=full_reduce, y=circuit) + with_assert and mode == STEINER_MODE and self.assertGates(circuit, architecture) with_assert and full_reduce and self.assertCircuitEquivalentNdArr(circuit, array) return circuit, matrix, rank @@ -174,7 +209,7 @@ def test_pso_optimization(self): with self.subTest(i=(p, mode, permute_input, permute_output)): circuits, perms, score = sequential_gauss([Mat2(np.copy(m)) for m in matrices], mode=mode, architecture=self.arch, full_reduce=True, n_steps=5, swarm_size=10, population_size=5, n_iterations=5, input_perm=permute_input, output_perm=permute_output) # It doesn't need to find an optimized solution, it only needs to do a non-trivial run - print(mode, score) + #print(mode, score) if not permute_input: self.assertListEqual(perms[0].tolist(), [i for i in range(self.n_qubits)]) if not permute_output: