HLS with coupling map

qiskit/circuit/library/generalized_gates/linear_function.py

@@ -171,6 +171,18 @@ def permutation_pattern(self):
         locs = np.where(linear == 1)
         return locs[1]
 
+    def permute(self, perm):
+        """Returns a linear function obtained by permuting rows and columns of a given linear function
+        using the permutation pattern ``perm``.
+        """
+        mat = self.linear
+        nq = mat.shape[0]
+        permuted_mat = np.zeros(mat.shape, dtype=bool)
+        for i in range(nq):
+            for j in range(nq):
+                permuted_mat[i, j] = mat[perm[i], perm[j]]
+        return LinearFunction(permuted_mat)
+
 
 def _linear_quantum_circuit_to_mat(qc: QuantumCircuit):
     """This creates a n x n matrix corresponding to the given linear quantum circuit."""

qiskit/synthesis/linear/linear_circuits_utils.py

@@ -13,9 +13,9 @@
 """Utility functions for handling linear reversible circuits."""
 
 import copy
-from typing import Callable
+from typing import Callable, List
 import numpy as np
-from qiskit import QuantumCircuit
+from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.exceptions import QiskitError
 from qiskit.circuit.exceptions import CircuitError
 from . import calc_inverse_matrix, check_invertible_binary_matrix
@@ -61,10 +61,25 @@ def optimize_cx_4_options(function: Callable, mat: np.ndarray, optimize_count: b
     if not check_invertible_binary_matrix(mat):
         raise QiskitError("The matrix is not invertible.")
 
+    circuits = _cx_circuits_4_options(function, mat)
+    best_qc = _choose_best_circuit(circuits, optimize_count)
+    return best_qc
+
+
+def _cx_circuits_4_options(function: Callable, mat: np.ndarray) -> List[QuantumCircuit]:
+    """Construct different circuits implementing a binary invertible matrix M,
+    by considering all four options: M,M^(-1),M^T,M^(-1)^T.
+
+    Args:
+        function: the synthesis function.
+        mat: a binary invertible matrix.
+
+    Returns:
+        List[QuantumCircuit]: constructed circuits.
+    """
+    circuits = []
     qc = function(mat)
-    best_qc = qc
-    best_depth = qc.depth()
-    best_count = qc.count_ops()["cx"]
+    circuits.append(qc)
 
     for i in range(1, 4):
         mat_cpy = copy.deepcopy(mat)
@@ -73,30 +88,67 @@ def optimize_cx_4_options(function: Callable, mat: np.ndarray, optimize_count: b
             mat_cpy = calc_inverse_matrix(mat_cpy)
             qc = function(mat_cpy)
             qc = qc.inverse()
+            circuits.append(qc)
         elif i == 2:
             mat_cpy = np.transpose(mat_cpy)
             qc = function(mat_cpy)
             qc = transpose_cx_circ(qc)
+            circuits.append(qc)
         elif i == 3:
             mat_cpy = calc_inverse_matrix(np.transpose(mat_cpy))
             qc = function(mat_cpy)
             qc = transpose_cx_circ(qc)
             qc = qc.inverse()
+            circuits.append(qc)
+
+    return circuits
 
-        new_depth = qc.depth()
-        new_count = qc.count_ops()["cx"]
-        # Prioritize count, and if it has the same count, then also consider depth
-        better_count = (optimize_count and best_count > new_count) or (
-            not optimize_count and best_depth == new_depth and best_count > new_count
-        )
-        # Prioritize depth, and if it has the same depth, then also consider count
-        better_depth = (not optimize_count and best_depth > new_depth) or (
-            optimize_count and best_count == new_count and best_depth > new_depth
-        )
-
-        if better_count or better_depth:
-            best_count = new_count
-            best_depth = new_depth
-            best_qc = qc
 
+def _choose_best_circuit(
+    circuits: List[QuantumCircuit], optimize_count: bool = True
+) -> QuantumCircuit:
+    """Returns the best quantum circuit either in terms of gate count or depth.
+
+    Args:
+        circuits: a list of quantum circuits
+        optimize_count: True if the number of CX gates is optimized, False if the depth is optimized.
+
+    Returns:
+        QuantumCircuit: the best quantum circuit out of the given circuits.
+    """
+    best_qc = circuits[0]
+    for circuit in circuits[1:]:
+        if _compare_circuits(best_qc, circuit, optimize_count=optimize_count):
+            best_qc = circuit
     return best_qc
+
+
+def _compare_circuits(
+    qc1: QuantumCircuit, qc2: QuantumCircuit, optimize_count: bool = True
+) -> bool:
+    """Compares two quantum circuits either in terms of gate count or depth.
+
+     Args:
+        qc1: the first quantum circuit
+        qc2: the second quantum circuit
+        optimize_count: True if the number of CX gates is optimized, False if the depth is optimized.
+
+    Returns:
+        bool: ``False`` means that the first quantum circuit is "better", ``True`` means the second.
+    """
+    # TODO: this is not fully correct if there are SWAP gates
+    count1 = qc1.size()
+    depth1 = qc1.depth()
+    count2 = qc2.size()
+    depth2 = qc2.depth()
+
+    # Prioritize count, and if it has the same count, then also consider depth
+    count2_is_better = (optimize_count and count1 > count2) or (
+        not optimize_count and depth1 == depth2 and count1 > count2
+    )
+    # Prioritize depth, and if it has the same depth, then also consider count
+    depth2_is_better = (not optimize_count and depth1 > depth2) or (
+        optimize_count and count1 == count2 and depth1 > depth2
+    )
+
+    return count2_is_better or depth2_is_better

qiskit/transpiler/passes/synthesis/high_level_synthesis.py

@@ -13,13 +13,16 @@
 
 """Synthesize higher-level objects."""
 
+from typing import Union, List
 
 from qiskit.converters import circuit_to_dag
 from qiskit.transpiler.basepasses import TransformationPass
 from qiskit.dagcircuit.dagcircuit import DAGCircuit
 from qiskit.transpiler.exceptions import TranspilerError
-from qiskit.synthesis import synth_clifford_full
-from qiskit.synthesis.linear import synth_cnot_count_full_pmh
+from qiskit.transpiler.coupling import CouplingMap
+from qiskit.synthesis.clifford import synth_clifford_full
+from qiskit.synthesis.linear import synth_cnot_count_full_pmh, synth_cnot_depth_line_kms
+from qiskit.synthesis.linear.linear_circuits_utils import optimize_cx_4_options, _compare_circuits
 from .plugin import HighLevelSynthesisPluginManager, HighLevelSynthesisPlugin
 
 
@@ -89,9 +92,20 @@ class HighLevelSynthesis(TransformationPass):
     ``default`` methods for all other high-level objects, including ``op_a``-objects.
     """
 
-    def __init__(self, hls_config=None):
+    def __init__(self, coupling_map: CouplingMap = None, hls_config=None):
+        """
+        HighLevelSynthesis initializer.
+
+        Args:
+            coupling_map (CouplingMap): the coupling map of the backend
+                in case synthesis is done on a physical circuit.
+            hls_config (HLSConfig): the high-level-synthesis config file
+            specifying synthesis methods and parameters.
+        """
         super().__init__()
 
+        self._coupling_map = coupling_map
+
         if hls_config is not None:
             self.hls_config = hls_config
         else:
@@ -110,9 +124,10 @@ def run(self, dag: DAGCircuit) -> DAGCircuit:
         Raises:
             TranspilerError: when the specified synthesis method is not available.
         """
-
         hls_plugin_manager = HighLevelSynthesisPluginManager()
 
+        dag_bit_indices = {bit: i for i, bit in enumerate(dag.qubits)}
+
         for node in dag.op_nodes():
 
             if node.name in self.hls_config.methods.keys():
@@ -141,14 +156,24 @@ def run(self, dag: DAGCircuit) -> DAGCircuit:
 
                 plugin_method = hls_plugin_manager.method(node.name, plugin_name)
 
-                # ToDo: similarly to UnitarySynthesis, we should pass additional parameters
-                #       e.g. coupling_map to the synthesis algorithm.
+                if self._coupling_map:
+                    plugin_args["coupling_map"] = self._coupling_map
+                    plugin_args["qubits"] = [dag_bit_indices[x] for x in node.qargs]
+
                 decomposition = plugin_method.run(node.op, **plugin_args)
 
-                # The synthesis methods that are not suited for the given higher-level-object
-                # will return None, in which case the next method in the list will be used.
+                # The above synthesis method may return:
+                # - None, if the synthesis algorithm is not suited for the given higher-level-object.
+                # - decomposition when the order of qubits is not important
+                # - a tuple (decomposition, wires) when the node's qubits need to be reordered
                 if decomposition is not None:
-                    dag.substitute_node_with_dag(node, circuit_to_dag(decomposition))
+                    if isinstance(decomposition, tuple):
+                        decomposition_dag = circuit_to_dag(decomposition[0])
+                        wires = [decomposition_dag.wires[i] for i in decomposition[1]]
+                        dag.substitute_node_with_dag(node, decomposition_dag, wires=wires)
+                    else:
+                        dag.substitute_node_with_dag(node, circuit_to_dag(decomposition))
+
                     break
 
         return dag
@@ -168,5 +193,160 @@ class DefaultSynthesisLinearFunction(HighLevelSynthesisPlugin):
 
     def run(self, high_level_object, **options):
         """Run synthesis for the given LinearFunction."""
-        decomposition = synth_cnot_count_full_pmh(high_level_object.linear)
+        # For now, use PMH algorithm by default
+        decomposition = PMHSynthesisLinearFunction().run(high_level_object, **options)
         return decomposition
+
+
+class KMSSynthesisLinearFunction(HighLevelSynthesisPlugin):
+    """Linear function synthesis plugin based on the Kutin-Moulton-Smithline method."""
+
+    def run(self, high_level_object, **options):
+        """Run synthesis for the given LinearFunction."""
+
+        # options supported by this plugin
+        coupling_map = options.get("coupling_map", None)
+        qubits = options.get("qubits", None)
+        consider_all_mats = options.get("all_mats", 0)
+        max_paths = options.get("max_paths", 1)
+        consider_original_circuit = options.get("orig_circuit", 1)
+        optimize_count = options.get("opt_count", 1)
+
+        # At the end, if not none, represents the best decomposition adhering to LNN architecture.
+        best_decomposition = None
+
+        # At the end, if not none, represents the path of qubits through the coupling map
+        # over which the LNN synthesis is applied.
+        best_path = None
+
+        if consider_original_circuit:
+            best_decomposition = high_level_object.original_circuit
+
+        if not coupling_map:
+            if not consider_all_mats:
+                decomposition = synth_cnot_depth_line_kms(high_level_object.linear)
+            else:
+                decomposition = optimize_cx_4_options(
+                    synth_cnot_depth_line_kms,
+                    high_level_object.linear,
+                    optimize_count=optimize_count,
+                )
+
+            if not best_decomposition or _compare_circuits(
+                best_decomposition, decomposition, optimize_count=optimize_count
+            ):
+                best_decomposition = decomposition
+
+        else:
+            # Consider the coupling map over the qubits on which the linear function is applied.
+            reduced_map = coupling_map.reduce(qubits)
+
+            # Find one or more paths through the coupling map (when such exist).
+            considered_paths = _hamiltonian_paths(reduced_map, max_paths)
+
+            for path in considered_paths:
+                permuted_linear_function = high_level_object.permute(path)
+
+                if not consider_all_mats:
+                    decomposition = synth_cnot_depth_line_kms(permuted_linear_function.linear)
+                else:
+                    decomposition = optimize_cx_4_options(
+                        synth_cnot_depth_line_kms,
+                        permuted_linear_function.linear,
+                        optimize_count=False,
+                    )
+
+                if not best_decomposition or _compare_circuits(
+                    best_decomposition, decomposition, optimize_count=False
+                ):
+                    best_decomposition = decomposition
+                    best_path = path
+
+        if best_path is None:
+            return best_decomposition
+
+        return best_decomposition, best_path
+
+
+class PMHSynthesisLinearFunction(HighLevelSynthesisPlugin):
+    """Linear function synthesis plugin based on the Patel-Markov-Hayes method."""
+
+    def run(self, high_level_object, **options):
+        """Run synthesis for the given LinearFunction."""
+
+        # options supported by this plugin
+        coupling_map = options.get("coupling_map", None)
+        consider_all_mats = options.get("all_mats", 0)
+        consider_original_circuit = options.get("orig_circuit", 1)
+        optimize_count = options.get("opt_count", 1)
+
+        # At the end, if not none, represents the best decomposition.
+        best_decomposition = None
+
+        if consider_original_circuit:
+            best_decomposition = high_level_object.original_circuit
+
+        # This synthesis method is not aware of the coupling map, so we cannot apply
+        # this method when the coupling map is not None.
+        # (Though, technically, we could check if the reduced coupling map is
+        # fully-connected).
+
+        if not coupling_map:
+            if not consider_all_mats:
+                decomposition = synth_cnot_count_full_pmh(high_level_object.linear)
+            else:
+                decomposition = optimize_cx_4_options(
+                    synth_cnot_count_full_pmh,
+                    high_level_object.linear,
+                    optimize_count=optimize_count,
+                )
+
+            if not best_decomposition or _compare_circuits(
+                best_decomposition, decomposition, optimize_count=optimize_count
+            ):
+                best_decomposition = decomposition
+
+        return best_decomposition
+
+
+def _hamiltonian_paths(
+    coupling_map: CouplingMap, cutoff: Union[None | int] = None
+) -> List[List[int]]:
+    """Returns a list of all Hamiltonian paths in ``coupling_map`` (stopping the enumeration when
+    the number of already discovered paths exceeds the ``cutoff`` value, when specified).
+    In particular, returns an empty list if there are no Hamiltonian paths.
+    """
+
+    # This is a temporary function, the plan is to move it to rustworkx
+
+    def should_stop():
+        return cutoff is not None and len(all_paths) >= cutoff
+
+    def _recurse(current_node):
+        current_path.append(current_node)
+        if len(current_path) == coupling_map.size():
+            # Discovered a new Hamiltonian path
+            all_paths.append(current_path.copy())
+
+        if should_stop():
+            return
+
+        unvisited_neighbors = [
+            node for node in coupling_map.neighbors(current_node) if node not in current_path
+        ]
+        for node in unvisited_neighbors:
+            _recurse(node)
+            if should_stop():
+                return
+
+        current_path.pop()
+
+    all_paths = []
+    current_path = []
+    qubits = coupling_map.physical_qubits
+
+    for qubit in qubits:
+        _recurse(qubit)
+        if should_stop():
+            break
+    return all_paths