1000x faster swap mapping and Cython build chain (#1789)

* Add Cython files.

* by error

* Modify swap mapper

* cleanup setup.py

* move cythonize to main setup.py

* fix style and add cython to requirements

* have travis build cython before testing

* disable lint checks that fail due to compiled code

* add cython to requirements (the txt one)

* revert to NumPy RandomState RNG for deterministic random generation

* remove unused import

* update tests

* make Cython build more portable

* move cython to dev requirements

* remove unused code

* turn on warnings, and fix old docstring

* fix import errors

* add cython build to appveyor

* add cython to appveyor

* revert aer changes

* fix "ground truth"

* this should not be here

* this should not be here

* this should not be here

* this should not be here

* this should not be here

* this should not be here

* I have no idea why they are saying they are modified

* Add Cython to setup_requires for setup.py

The setup_requires parameter to the setup() function is used to
specify packages that need to be present in order for the setup script
to run. [1] With the introduction of Cython to the setup.py this is now
required to be installed for setup, so this commit adds Cython to
setup_requires to indicate this.

[1] https://setuptools.readthedocs.io/en/latest/setuptools.html

* add to changelog

* Update Cython version in setup_requires

.travis.yml

@@ -32,6 +32,7 @@ stage_generic: &stage_generic
     - pip install qiskit-ibmq-provider
   script:
     # Compile the executables and run the tests.
+    - python setup.py build_ext --inplace
     - make test_ci
   after_failure:
     - python tools/report_ci_failure.py

CHANGELOG.rst

@@ -23,7 +23,7 @@ The format is based on `Keep a Changelog`_.
 Added
 -----
 
-
+- Core StochasticSwap routine implimented in Cython (#1789).
 - New EnlargeWithAncilla pass for adding ancilla qubits after a Layout
   selection pass (#1603).
 - New Unroll2Q pass for unrolling gates down to just 1q or 2q gates (#1614).

appveyor.yml

@@ -34,9 +34,11 @@ install:
     - pip.exe install vcrpy
     - pip.exe install jupyter
     - pip.exe install ipywidgets
+    - pip.exe install cython
     - pip.exe install stestr
     - pip.exe install qiskit-aer
     - pip.exe install qiskit-ibmq-provider
+    - python setup.py build_ext --inplace
 # TODO: uncomment this when testing-cabal/subunit#33 is fixed
 #   - pip.exe install junitxml
 

qiskit/transpiler/passes/mapping/cython/__init__.py

@@ -0,0 +1,8 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2018, IBM.
+#
+# This source code is licensed under the Apache License, Version 2.0 found in
+# the LICENSE.txt file in the root directory of this source tree.
+
+"""Module containing transpiler Cython code."""

qiskit/transpiler/passes/mapping/cython/stochastic_swap/__init__.py

@@ -0,0 +1,8 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2018, IBM.
+#
+# This source code is licensed under the Apache License, Version 2.0 found in
+# the LICENSE.txt file in the root directory of this source tree.
+
+"""Module containing Cython code for StochasticSwap mapper."""

qiskit/transpiler/passes/mapping/cython/stochastic_swap/_swap_trial.pyx

@@ -0,0 +1,186 @@
+# -*- coding: utf-8 -*-
+#!python
+#cython: language_level = 3
+#distutils: language = c++
+
+# Copyright 2018, IBM.
+#
+# This source code is licensed under the Apache License, Version 2.0 found in
+# the LICENSE.txt file in the root directory of this source tree.
+
+cimport cython
+from libcpp.set cimport set as cset
+from .utils cimport NLayout, EdgeCollection
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef double compute_cost(double[:, ::1] dist, unsigned int * logic_to_phys,
+                          int[::1] gates, unsigned int num_gates) nogil:
+    """ Computes the cost (distance) of a logical to physical mapping.
+    
+    Args:
+        dist (ndarray): An array of doubles that specifies the distance.
+        logic_to_phys (int *): Pointer to logical to physical array.
+        gates (ndarray): Array of ints giving gates in layer.
+        num_gates (int): The number of gates (length of gates//2).
+    
+    Returns:
+        double: The distance calculated.
+    """
+    cdef unsigned int ii, jj, kk
+    cdef double cost = 0.0
+    for kk in range(num_gates):
+        ii = logic_to_phys[gates[2*kk]]
+        jj = logic_to_phys[gates[2*kk+1]]
+        cost += dist[ii,jj]
+    return cost
+
+@cython.nonecheck(False)
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef compute_random_scaling(double[:, ::1] scale, double[:, ::1] cdist2,
+                            double * rand, int num_qubits):
+    """ Computes the symmetric random scaling (perturbation) matrix, 
+    and places the values in the 'scale' array.
+
+    Args:
+        scale (ndarray): An array of doubles where the values are to be stored.
+        cdist2 (ndarray): Array representing the coupling map distance squared.
+        rand (double *): Array of rands of length num_qubits*(num_qubits+1)//2.
+        num_qubits (int): Number of physical qubits.
+    """
+    cdef size_t ii, jj, idx=0
+    for ii in range(num_qubits):
+        for jj in range(ii):
+            scale[ii,jj] = rand[idx]*cdist2[ii,jj]
+            scale[jj,ii] = scale[ii,jj]
+            idx += 1
+
+
+@cython.nonecheck(False)
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def swap_trial(int num_qubits, NLayout int_layout, int[::1] int_qubit_subset,
+               int[::1] gates, double[:, ::1] cdist2, double[:, ::1] cdist, 
+               int[::1] edges, double[:, ::1] scale, object rng):
+    """ A single iteration of the tchastic swap mapping routine.
+
+    Args:
+        num_qubits (int): The number of physical qubits.
+        int_layout (NLayout): The numeric (integer) representation of 
+                              the initial_layout.
+        int_qubit_subset (ndarray): Int ndarray listing qubits in set.
+        gates (ndarray): Int array with integers giving qubits on which
+                         two-qubits gates act on.
+        cdist2 (ndarray): Array of doubles that gives the square of the 
+                          distance graph.
+        cdist (ndarray): Array of doubles that gives the distance graph.
+        edges (ndarray): Int array of edges in coupling map.
+        scale (ndarray): A double array that holds the perturbed cdist2 array.
+        rng (RandomState): An instance of the NumPy RandomState.
+
+    Returns:
+        double: Best distance achieved in this trial.
+        EdgeCollection: Collection of optimal edges found.
+        NLayout: The optimal layout found.
+        int: The number of depth steps required in mapping.
+    """
+    cdef EdgeCollection opt_edges = EdgeCollection()
+    cdef NLayout optimal_layout, new_layout, trial_layout = int_layout.copy()
+    
+    cdef unsigned int num_gates = gates.shape[0]//2
+    cdef unsigned int num_edges = edges.shape[0]//2
+    
+    cdef unsigned int need_copy, cost_reduced
+    cdef unsigned int depth_step = 1
+    cdef unsigned int depth_max = 2 * num_qubits + 1
+    cdef double min_cost, new_cost, dist
+    
+    cdef unsigned int start_edge, end_edge, start_qubit, end_qubit
+    cdef unsigned int optimal_start, optimal_end, optimal_start_qubit, optimal_end_qubit
+    
+    cdef size_t idx
+    
+    # Compute randomized distance
+    cdef double[::1] rand = 1.0 + rng.normal(0.0, 1.0/num_qubits,
+                                             size=num_qubits*(num_qubits+1)//2)
+    
+    compute_random_scaling(scale, cdist2, &rand[0], num_qubits)
+    
+    # Convert int qubit array to c++ set
+    cdef cset[unsigned int] qubit_set
+    cdef cset[unsigned int] input_qubit_set
+    
+    for idx in range(int_qubit_subset.shape[0]):
+        input_qubit_set.insert(int_qubit_subset[idx])
+    
+    # Loop over depths from 1 up to a maximum depth
+    while depth_step < depth_max:
+        qubit_set = input_qubit_set
+        # While there are still qubits available
+        while not qubit_set.empty():
+            # Compute the objective function
+            min_cost = compute_cost(scale, trial_layout.logic_to_phys,
+                                   gates, num_gates)
+            # Try to decrease objective function
+            cost_reduced = 0
+
+            # Loop over edges of coupling graph
+            need_copy = 1
+            for idx in range(num_edges):
+                start_edge = edges[2*idx]
+                end_edge = edges[2*idx+1]
+                start_qubit = trial_layout.phys_to_logic[start_edge]
+                end_qubit =  trial_layout.phys_to_logic[end_edge]
+                # Are the qubits available?
+                if  qubit_set.count(start_qubit) and qubit_set.count(end_qubit):
+                    # Try this edge to reduce the cost
+                    if need_copy:
+                        new_layout = trial_layout.copy()
+                        need_copy = 0
+                    new_layout.swap(start_edge, end_edge)
+                    # Compute the objective function
+                    new_cost = compute_cost(scale, new_layout.logic_to_phys,
+                                   gates, num_gates)
+                    # Record progress if we succceed
+                    if new_cost < min_cost:
+                        cost_reduced = True
+                        min_cost = new_cost
+                        optimal_layout = new_layout
+                        optimal_start = start_edge
+                        optimal_end = end_edge
+                        optimal_start_qubit = start_qubit
+                        optimal_end_qubit = end_qubit
+                        need_copy = 1
+                    else:
+                        new_layout.swap(start_edge, end_edge)
+
+            # After going over all edges
+            # Were there any good swap choices?
+            if cost_reduced:
+                qubit_set.erase(optimal_start_qubit)
+                qubit_set.erase(optimal_end_qubit)
+                trial_layout = optimal_layout
+                opt_edges.add(optimal_start, optimal_end)
+            else:
+                break
+
+        # We have either run out of swap pairs to try or
+        # failed to improve the cost.
+
+        # Compute the coupling graph distance
+        dist = compute_cost(cdist, trial_layout.logic_to_phys,
+                                   gates, num_gates)
+        # If all gates can be applied now, we are finished.
+        # Otherwise we need to consider a deeper swap circuit
+        if dist == num_gates:
+            break
+
+        # Increment the depth
+        depth_step += 1
+
+    # Either we have succeeded at some depth d < dmax or failed
+    dist = compute_cost(cdist, trial_layout.logic_to_phys,
+                                   gates, num_gates)
+    
+    return dist, opt_edges, trial_layout, depth_step

qiskit/transpiler/passes/mapping/cython/stochastic_swap/utils.pxd

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+#!python
+#cython: language_level = 3, cdivision = True, nonecheck = False
+#distutils: language = c++
+
+# Copyright 2018, IBM.
+#
+# This source code is licensed under the Apache License, Version 2.0 found in
+# the LICENSE.txt file in the root directory of this source tree.
+from libcpp.vector cimport vector
+
+# Numeric layout --------------------------------------------------------------
+cdef class NLayout:
+    cdef:
+        unsigned int l2p_len
+        unsigned int p2l_len
+        unsigned int * logic_to_phys
+        unsigned int * phys_to_logic
+
+    # Methods
+    cdef NLayout copy(self)
+    cdef void swap(self, unsigned int idx1, unsigned int idx2)
+    cpdef object to_layout(self, object dag)
+
+
+cpdef NLayout nlayout_from_layout(object layout,
+                                  object dag, 
+                                  unsigned int physical_qubits)
+
+
+# Edge collection -------------------------------------------------------------
+cdef class EdgeCollection:
+    cdef vector[unsigned int] _edges
+
+    cpdef void add(self, unsigned int edge_start, unsigned int edge_end)

qiskit/transpiler/passes/mapping/cython/stochastic_swap/utils.pyx

@@ -0,0 +1,188 @@
+# -*- coding: utf-8 -*-
+#!python
+#cython: language_level = 3
+#distutils: language = c++
+
+# Copyright 2018, IBM.
+#
+# This source code is licensed under the Apache License, Version 2.0 found in
+# the LICENSE.txt file in the root directory of this source tree.
+
+cimport cython
+import numpy as np
+from libc.stdlib cimport calloc, free
+from libcpp.vector cimport vector
+
+from qiskit.mapper.layout import Layout
+
+
+cdef class EdgeCollection:
+    """ A simple contain that contains a C++ vector
+    representing edges in the coupling map that are
+    found to be optimal by the swap mapper.  This allows
+    us to keep the vector alive.
+    """
+    cpdef void add(self, unsigned int edge_start, unsigned int edge_end):
+        """ Add two edges, in order, to the collection.
+
+        Args:
+            edge_start (int): The beginning edge.
+            edge_end (int): The end of the edge.
+        """
+        self._edges.push_back(edge_start)
+        self._edges.push_back(edge_end)
+    
+    @property
+    def size(self):
+        """ The size of the edge collection.
+        Returns:
+            int: Size of the edge collection.
+        """
+        return self._edges.size()
+    @cython.boundscheck(False)
+    def edges(self):
+        """ Returns the vector of edges as a NumPy arrau.
+        Returns:
+            ndarray: Int array of edges.
+        """
+        cdef size_t kk
+        out = np.zeros(self._edges.size(), dtype=np.uint32)
+        for kk in range(self._edges.size()):
+            out[kk] = self._edges[kk]
+        return out
+
+
+cdef class NLayout:
+    """ A Numeric representation of a Qiskit Layout object.
+    Here all qubit layouts are stored as int arrays.
+    """
+    def __cinit__(self, unsigned int num_logical,
+                  unsigned int num_physical):
+        """ Init object.
+        Args:
+            num_logical (int): Number of logical qubits.
+            num_physical (int): Number of physical qubits.
+        """
+        self.l2p_len = num_logical
+        self.p2l_len = num_physical
+        self.logic_to_phys = <unsigned int *>calloc(num_logical,
+                                                    sizeof(unsigned int))
+        self.phys_to_logic = <unsigned int *>calloc(num_physical,
+                                                    sizeof(unsigned int))
+    
+    def __dealloc__(self):
+        """ Clears the pointers when finished.
+        """
+        if self.logic_to_phys is not NULL:
+            free(self.logic_to_phys)
+            self.logic_to_phys = NULL
+        if self.phys_to_logic is not NULL:
+            free(self.phys_to_logic)
+            self.phys_to_logic = NULL
+            
+    @property
+    def logic_to_phys(self):
+        """ The array mapping logical to physical qubits.
+        Returns:
+            ndarray: Int array of logical to physical mappings.
+        """
+        cdef size_t kk
+        out = np.zeros(self.l2p_len, dtype=np.int32)
+        for kk in range(self.l2p_len):
+            out[kk] = self.logic_to_phys[kk]
+        return out
+    
+    @property
+    def phys_to_logic(self):
+        """ The array mapping physical to logical qubits.
+        Returns:
+            ndarray: Int array of physical to logical mappings.
+        """
+        cdef size_t kk
+        out = np.zeros(self.p2l_len, dtype=np.int32)
+        for kk in range(self.p2l_len):
+            out[kk] = self.phys_to_logic[kk]
+        return out
+    
+    @cython.boundscheck(False)
+    cdef NLayout copy(self):
+        """ Returns a copy of the layout.
+
+        Returns:
+            NLayout: A copy of the layout.
+        """
+        cdef NLayout out = NLayout(self.l2p_len, self.p2l_len)
+        cdef size_t kk
+        for kk in range(self.l2p_len):
+            out.logic_to_phys[kk] = self.logic_to_phys[kk]
+        for kk in range(self.p2l_len):
+            out.phys_to_logic[kk] = self.phys_to_logic[kk]
+        return out
+            
+    @cython.boundscheck(False)
+    cdef void swap(self, unsigned int idx1, unsigned int idx2):
+        """ Swaps two indices in the Layout
+
+        Args:
+            idx1 (int): Index 1.
+            idx2 (int): Index 2.
+        """
+        cdef unsigned int temp1, temp2
+        temp1 = self.phys_to_logic[idx1]
+        temp2 = self.phys_to_logic[idx2]
+        self.phys_to_logic[idx1] = temp2
+        self.phys_to_logic[idx2] = temp1
+        self.logic_to_phys[self.phys_to_logic[idx1]] = idx1
+        self.logic_to_phys[self.phys_to_logic[idx2]] = idx2
+        
+    @cython.boundscheck(False)
+    cpdef object to_layout(self, object qregs):
+        """ Converts numeric layout back to Qiskit Layout object.
+
+        Args:
+            qregs (OrderedDict): An ordered dict of (QuantumRegister, int)
+                                 tuples.
+        
+        Returns:
+            Layout: The corresponding Qiskit Layout object.
+        """
+        out = Layout()
+        cdef unsigned int main_idx = 0
+        cdef size_t idx
+        for qreg in qregs.values():
+            for idx in range(qreg.size):
+                out[(qreg, idx)] = self.logic_to_phys[main_idx]
+                main_idx += 1
+        return out
+    
+    
+cpdef NLayout nlayout_from_layout(object layout, object qregs, 
+                                  unsigned int physical_qubits):
+    """ Converts Qiskit Layout object to numerical NLayout.
+
+    Args:
+        layout (Layout): A Qiskit Layout instance.
+        qregs (OrderedDict): An ordered dict of (QuantumRegister, int)
+                                tuples.
+        physical_qubits (int): Number of physical qubits.
+    Returns:
+        NLayout: The corresponding numerical layout.
+    """
+    cdef size_t ind
+    cdef list sizes = [qr.size for qr in qregs.values()]
+    cdef int[::1] reg_idx = np.cumsum([0]+sizes, dtype=np.int32)
+    cdef unsigned int logical_qubits = sum(sizes)
+
+    cdef dict regint = {}
+    for ind, qreg in enumerate(qregs.values()):
+        regint[qreg] = ind
+
+    cdef NLayout out = NLayout(logical_qubits, physical_qubits)
+    cdef object key, val
+    cdef dict merged_dict = {**layout._p2v, **layout._v2p}
+    for key, val in merged_dict.items():
+        if isinstance(key, tuple):
+            out.logic_to_phys[reg_idx[regint[key[0]]]+key[1]] = val
+        else:
+            out.phys_to_logic[key] = reg_idx[regint[val[0]]]+val[1]
+    return out

qiskit/transpiler/passes/mapping/stochastic_swap.py

@@ -20,7 +20,10 @@ from qiskit.dagcircuit import DAGCircuit
 from qiskit.extensions.standard import SwapGate
 from qiskit.mapper import Layout
 from .barrier_before_final_measurements import BarrierBeforeFinalMeasurements
-
+# pylint: disable=no-name-in-module, import-error
+from .cython.stochastic_swap.utils import nlayout_from_layout
+# pylint: disable=no-name-in-module, import-error
+from .cython.stochastic_swap._swap_trial import swap_trial
 logger = getLogger(__name__)
 
 
@@ -72,6 +75,8 @@ class StochasticSwap(TransformationPass):
         self.input_layout = None
         self.trials = trials
         self.seed = seed
+        self.qregs = None
+        self.rng = None
         self.requires.append(BarrierBeforeFinalMeasurements())
 
     def run(self, dag):
@@ -104,12 +109,15 @@ class StochasticSwap(TransformationPass):
 
         self.input_layout = self.initial_layout.copy()
 
-        new_dag = self._mapper(dag, self.coupling_map, trials=self.trials, seed=self.seed)
+        self.qregs = dag.qregs
+        self.rng = np.random.RandomState(self.seed)
+
+        new_dag = self._mapper(dag, self.coupling_map, trials=self.trials)
         # self.property_set["layout"] = self.initial_layout
         return new_dag
 
     def _layer_permutation(self, layer_partition, layout, qubit_subset,
-                           coupling, trials, seed=None):
+                           coupling, trials):
         """Find a swap circuit that implements a permutation for this layer.
 
         The goal is to swap qubits such that qubits in the same two-qubit gates
@@ -129,11 +137,9 @@ class StochasticSwap(TransformationPass):
             This coupling map should be one that was provided to the
             stochastic mapper.
         trials (int): Number of attempts the randomized algorithm makes.
-        seed (int): Optional seed for the random number generator. If it is
-            None we do not reseed.
 
         Returns:
-             Tuple: success_flag, best_circuit, best_depth, best_layout, trivial_flag
+            Tuple: success_flag, best_circuit, best_depth, best_layout, trivial_flag
 
         If success_flag is True, then best_circuit contains a DAGCircuit with
         the swap circuit, best_depth contains the depth of the swap circuit,
@@ -143,168 +149,11 @@ class StochasticSwap(TransformationPass):
 
         Raises:
             TranspilerError: if anything went wrong.
-        """
-        if seed is not None:
-            np.random.seed(seed)
-
-        logger.debug("layer_permutation: layer_partition = %s",
-                     pformat(layer_partition))
-        logger.debug("layer_permutation: layout = %s",
-                     pformat(layout.get_virtual_bits()))
-        logger.debug("layer_permutation: qubit_subset = %s",
-                     pformat(qubit_subset))
-        logger.debug("layer_permutation: trials = %s", trials)
-
-        gates = []  # list of lists of tuples [[(register, index), ...], ...]
-        for gate_args in layer_partition:
-            if len(gate_args) > 2:
-                raise TranspilerError("Layer contains > 2-qubit gates")
-            elif len(gate_args) == 2:
-                gates.append(tuple(gate_args))
-        logger.debug("layer_permutation: gates = %s", pformat(gates))
-
-        # Can we already apply the gates? If so, there is no work to do.
-        dist = sum([coupling.distance(layout[g[0]], layout[g[1]])
-                    for g in gates])
-        logger.debug("layer_permutation: distance = %s", dist)
-        if dist == len(gates):
-            logger.debug("layer_permutation: nothing to do")
-            circ = DAGCircuit()
-            for register in layout.get_virtual_bits().keys():
-                if register[0] not in circ.qregs.values():
-                    circ.add_qreg(register[0])
-            return True, circ, 0, layout, (not bool(gates))
-
-        # Begin loop over trials of randomized algorithm
-        num_qubits = len(layout)
-        best_depth = inf  # initialize best depth
-        best_circuit = None  # initialize best swap circuit
-        best_layout = None  # initialize best final layout
-
-        cdist2 = coupling._dist_matrix**2
-        # Scaling matrix
-        scale = np.zeros((num_qubits, num_qubits))
-        utri_idx = np.triu_indices(num_qubits)
-
-        for trial in range(trials):
-            logger.debug("layer_permutation: trial %s", trial)
-            trial_layout = layout.copy()
-            trial_circuit = DAGCircuit()  # SWAP circuit for this trial
-            for register in trial_layout.get_virtual_bits().keys():
-                if register[0] not in trial_circuit.qregs.values():
-                    trial_circuit.add_qreg(register[0])
-
-            # Compute randomized distance
-            data = 1 + np.random.normal(0, 1/num_qubits,
-                                        size=num_qubits*(num_qubits+1)//2)
-            scale[utri_idx] = data
-            xi = (scale+scale.T)*cdist2  # pylint: disable=invalid-name
-
-            slice_circuit = DAGCircuit()  # circuit for this swap slice
-            for register in trial_layout.get_virtual_bits().keys():
-                if register[0] not in slice_circuit.qregs.values():
-                    slice_circuit.add_qreg(register[0])
-
-            # Loop over depths from 1 up to a maximum depth
-            depth_step = 1
-            depth_max = 2 * num_qubits + 1
-            while depth_step < depth_max:
-                qubit_set = set(qubit_subset)
-                # While there are still qubits available
-                while qubit_set:
-                    # Compute the objective function
-                    min_cost = sum(xi[trial_layout[g[0]]][trial_layout[g[1]]] for g in gates)
-                    # Try to decrease objective function
-                    cost_reduced = False
-
-                    # Loop over edges of coupling graph
-                    need_copy = True
-                    for edge in coupling.get_edges():
-                        qubits = (trial_layout[edge[0]], trial_layout[edge[1]])
-                        # Are the qubits available?
-                        if qubits[0] in qubit_set and qubits[1] in qubit_set:
-                            # Try this edge to reduce the cost
-                            if need_copy:
-                                new_layout = trial_layout.copy()
-                                need_copy = False
-                            new_layout.swap(edge[0], edge[1])
-
-                            # Compute the objective function
-                            new_cost = sum(xi[new_layout[g[0]]][new_layout[g[1]]] for g in gates)
-                            # Record progress if we succceed
-                            if new_cost < min_cost:
-                                logger.debug("layer_permutation: min_cost "
-                                             "improved to %s", min_cost)
-                                cost_reduced = True
-                                min_cost = new_cost
-                                optimal_layout = new_layout
-                                optimal_edge = (self.initial_layout[edge[0]],
-                                                self.initial_layout[edge[1]])
-                                optimal_qubits = qubits
-                                need_copy = True
-                            else:
-                                new_layout.swap(edge[0], edge[1])
-
-                    # Were there any good swap choices?
-                    if cost_reduced:
-                        qubit_set.remove(optimal_qubits[0])
-                        qubit_set.remove(optimal_qubits[1])
-                        trial_layout = optimal_layout
-                        slice_circuit.apply_operation_back(
-                            SwapGate(optimal_edge[0],
-                                     optimal_edge[1]))
-                        logger.debug("layer_permutation: swap the pair %s",
-                                     pformat(optimal_edge))
-                    else:
-                        break
-
-                # We have either run out of swap pairs to try or
-                # failed to improve the cost.
-
-                # Compute the coupling graph distance
-                dist = sum(coupling.distance(trial_layout[g[0]],
-                                             trial_layout[g[1]])
-                           for g in gates)
-                logger.debug("layer_permutation: new swap distance = %s", dist)
-                # If all gates can be applied now, we are finished.
-                # Otherwise we need to consider a deeper swap circuit
-                if dist == len(gates):
-                    logger.debug("layer_permutation: all gates can be "
-                                 "applied now in this layer")
-                    trial_circuit.extend_back(slice_circuit)
-                    break
-
-                # Increment the depth
-                depth_step += 1
-                logger.debug("layer_permutation: increment depth to %s", depth_step)
-
-            # Either we have succeeded at some depth d < dmax or failed
-            dist = sum(coupling.distance(trial_layout[g[0]],
-                                         trial_layout[g[1]])
-                       for g in gates)
-            logger.debug("layer_permutation: final distance for this trial = %s", dist)
-            if dist == len(gates):
-                if depth_step < best_depth:
-                    logger.debug("layer_permutation: got circuit with improved depth %s",
-                                 depth_step)
-                    best_circuit = trial_circuit
-                    best_layout = trial_layout
-                    best_depth = min(best_depth, depth_step)
-
-            # Break out of trial loop if we found a depth 1 circuit
-            # since we can't improve it further
-            if best_depth == 1:
-                break
-
-        # If we have no best circuit for this layer, all of the
-        # trials have failed
-        if best_circuit is None:
-            logger.debug("layer_permutation: failed!")
-            return False, None, None, None, False
-
-        # Otherwise, we return our result for this layer
-        logger.debug("layer_permutation: success!")
-        return True, best_circuit, best_depth, best_layout, False
+     """
+        return _layer_permutation(layer_partition, self.initial_layout,
+                                  layout, qubit_subset,
+                                  coupling, trials,
+                                  self.qregs, self.rng)
 
     def _layer_update(self, i, first_layer, best_layout, best_depth,
                       best_circuit, layer_list):
@@ -362,7 +211,7 @@ class StochasticSwap(TransformationPass):
         return dagcircuit_output
 
     def _mapper(self, circuit_graph, coupling_graph,
-                trials=20, seed=None):
+                trials=20):
         """Map a DAGCircuit onto a CouplingMap using swap gates.
 
         Use self.initial_layout for the initial layout.
@@ -371,7 +220,6 @@ class StochasticSwap(TransformationPass):
             circuit_graph (DAGCircuit): input DAG circuit
             coupling_graph (CouplingMap): coupling graph to map onto
             trials (int): number of trials.
-            seed (int): initial seed.
 
         Returns:
             DAGCircuit: object containing a circuit equivalent to
@@ -441,7 +289,7 @@ class StochasticSwap(TransformationPass):
             success_flag, best_circuit, best_depth, best_layout, trivial_flag \
                 = self._layer_permutation(layer["partition"], layout,
                                           qubit_subset, coupling_graph,
-                                          trials, seed)
+                                          trials)
             logger.debug("mapper: layer %d", i)
             logger.debug("mapper: success_flag=%s,best_depth=%s,trivial_flag=%s",
                          success_flag, str(best_depth), trivial_flag)
@@ -460,7 +308,7 @@ class StochasticSwap(TransformationPass):
                             serial_layer["partition"],
                             layout, qubit_subset,
                             coupling_graph,
-                            trials, seed)
+                            trials)
                     logger.debug("mapper: layer %d, sublayer %d", i, j)
                     logger.debug("mapper: success_flag=%s,best_depth=%s,"
                                  "trivial_flag=%s",
@@ -534,3 +382,171 @@ class StochasticSwap(TransformationPass):
                 dagcircuit_output.compose_back(layer["graph"], edge_map)
 
         return dagcircuit_output
+
+
+def _layer_permutation(layer_partition, initial_layout, layout, qubit_subset,
+                       coupling, trials, qregs, rng):
+    """Find a swap circuit that implements a permutation for this layer.
+
+    Args:
+        layer_partition (list): The layer_partition is a list of (qu)bit
+            lists and each qubit is a tuple (qreg, index).
+        initial_layout (Layout): The initial layout passed.
+        layout (Layout): The layout is a Layout object mapping virtual
+            qubits in the input circuit to physical qubits in the coupling
+            graph. It reflects the current positions of the data.
+        qubit_subset (list): The qubit_subset is the set of qubits in
+            the coupling graph that we have chosen to map into, as tuples
+            (Register, index).
+        coupling (CouplingMap): Directed graph representing a coupling map.
+            This coupling map should be one that was provided to the
+            stochastic mapper.
+        trials (int): Number of attempts the randomized algorithm makes.
+        qregs (OrderedDict): Ordered dict of registers from input DAG.
+        rng (RandomState): Random number generator.
+
+    Returns:
+        Tuple: success_flag, best_circuit, best_depth, best_layout, trivial_flag
+
+    Raises:
+        TranspilerError: if anything went wrong.
+     """
+    logger.debug("layer_permutation: layer_partition = %s",
+                 pformat(layer_partition))
+    logger.debug("layer_permutation: layout = %s",
+                 pformat(layout.get_virtual_bits()))
+    logger.debug("layer_permutation: qubit_subset = %s",
+                 pformat(qubit_subset))
+    logger.debug("layer_permutation: trials = %s", trials)
+
+    gates = []  # list of lists of tuples [[(register, index), ...], ...]
+    for gate_args in layer_partition:
+        if len(gate_args) > 2:
+            raise TranspilerError("Layer contains > 2-qubit gates")
+        elif len(gate_args) == 2:
+            gates.append(tuple(gate_args))
+    logger.debug("layer_permutation: gates = %s", pformat(gates))
+
+    # Can we already apply the gates? If so, there is no work to do.
+    dist = sum([coupling.distance(layout[g[0]], layout[g[1]])
+                for g in gates])
+    logger.debug("layer_permutation: distance = %s", dist)
+    if dist == len(gates):
+        logger.debug("layer_permutation: nothing to do")
+        circ = DAGCircuit()
+        for register in layout.get_virtual_bits().keys():
+            if register[0] not in circ.qregs.values():
+                circ.add_qreg(register[0])
+        return True, circ, 0, layout, (not bool(gates))
+
+    # Begin loop over trials of randomized algorithm
+    num_qubits = len(layout)
+    best_depth = inf  # initialize best depth
+    best_edges = None  # best edges found
+    best_circuit = None  # initialize best swap circuit
+    best_layout = None  # initialize best final layout
+
+    cdist2 = coupling._dist_matrix**2
+    # Scaling matrix
+    scale = np.zeros((num_qubits, num_qubits))
+
+    int_qubit_subset = regtuple_to_numeric(qubit_subset, qregs)
+    int_gates = gates_to_idx(gates, qregs)
+    int_layout = nlayout_from_layout(layout, qregs, coupling.size())
+
+    trial_circuit = DAGCircuit()  # SWAP circuit for this trial
+    for register in layout.get_virtual_bits().keys():
+        if register[0] not in trial_circuit.qregs.values():
+            trial_circuit.add_qreg(register[0])
+
+    slice_circuit = DAGCircuit()  # circuit for this swap slice
+    for register in layout.get_virtual_bits().keys():
+        if register[0] not in slice_circuit.qregs.values():
+            slice_circuit.add_qreg(register[0])
+    edges = np.asarray(coupling.get_edges(), dtype=np.int32).ravel()
+    cdist = coupling._dist_matrix
+    for trial in range(trials):
+        logger.debug("layer_permutation: trial %s", trial)
+        # This is one Trial --------------------------------------
+        dist, optim_edges, trial_layout, depth_step = swap_trial(num_qubits, int_layout,
+                                                                 int_qubit_subset,
+                                                                 int_gates, cdist2,
+                                                                 cdist, edges, scale,
+                                                                 rng)
+
+        logger.debug("layer_permutation: final distance for this trial = %s", dist)
+        if dist == len(gates) and depth_step < best_depth:
+            logger.debug("layer_permutation: got circuit with improved depth %s",
+                         depth_step)
+            best_edges = optim_edges
+            best_layout = trial_layout
+            best_depth = min(best_depth, depth_step)
+
+        # Break out of trial loop if we found a depth 1 circuit
+        # since we can't improve it further
+        if best_depth == 1:
+            break
+
+    # If we have no best circuit for this layer, all of the
+    # trials have failed
+    if best_layout is None:
+        logger.debug("layer_permutation: failed!")
+        return False, None, None, None, False
+
+    edgs = best_edges.edges()
+    for idx in range(best_edges.size//2):
+        slice_circuit.apply_operation_back(SwapGate(initial_layout[edgs[2*idx]],
+                                                    initial_layout[edgs[2*idx+1]]))
+    trial_circuit.extend_back(slice_circuit)
+    best_circuit = trial_circuit
+
+    # Otherwise, we return our result for this layer
+    logger.debug("layer_permutation: success!")
+    best_lay = best_layout.to_layout(qregs)
+    return True, best_circuit, best_depth, best_lay, False
+
+
+def regtuple_to_numeric(items, qregs):
+    """Takes (QuantumRegister, int) tuples and converts
+    them into an integer array.
+
+    Args:
+        items (list): List of tuples of (QuantumRegister, int)
+                      to convert.
+        qregs (dict): List of )QuantumRegister, int) tuples.
+    Returns:
+        ndarray: Array of integers.
+
+    """
+    sizes = [qr.size for qr in qregs.values()]
+    reg_idx = np.cumsum([0]+sizes)
+    regint = {}
+    for ind, qreg in enumerate(qregs.values()):
+        regint[qreg] = ind
+    out = np.zeros(len(items), dtype=np.int32)
+    for idx, val in enumerate(items):
+        out[idx] = reg_idx[regint[val[0]]]+val[1]
+    return out
+
+
+def gates_to_idx(gates, qregs):
+    """Converts gate tuples into a nested list of integers.
+
+    Args:
+        gates (list): List of (QuantumRegister, int) pairs
+                      representing gates.
+        qregs (dict): List of )QuantumRegister, int) tuples.
+
+    Returns:
+        list: Nested list of integers for gates.
+    """
+    sizes = [qr.size for qr in qregs.values()]
+    reg_idx = np.cumsum([0]+sizes)
+    regint = {}
+    for ind, qreg in enumerate(qregs.values()):
+        regint[qreg] = ind
+    out = np.zeros(2*len(gates), dtype=np.int32)
+    for idx, gate in enumerate(gates):
+        out[2*idx] = reg_idx[regint[gate[0][0]]]+gate[0][1]
+        out[2*idx+1] = reg_idx[regint[gate[1][0]]]+gate[1][1]
+    return out

requirements-dev.txt

@@ -11,3 +11,4 @@ stestr>=2.0.0
 vcrpy
 PyGithub
 wheel
+cython>=0.27.1

setup.py

@@ -5,10 +5,15 @@
 # This source code is licensed under the Apache License, Version 2.0 found in
 # the LICENSE.txt file in the root directory of this source tree.
 
-from setuptools import setup, find_packages
+"The Qiskit Terra setup file."
 
+import os
+import sys
+import distutils.sysconfig
+from setuptools import setup, find_packages, Extension
+from Cython.Build import cythonize
 
-requirements = [
+REQUIREMENTS = [
     "jsonschema>=2.6,<2.7",
     "marshmallow>=2.17.0,<3",
     "marshmallow_polyfield>=3.2,<4",
@@ -21,12 +26,49 @@ requirements = [
     "sympy>=1.3"
 ]
 
+# Add Cython extensions here
+CYTHON_EXTS = ['utils', '_swap_trial']
+CYTHON_MODULE = 'qiskit.transpiler.passes.mapping.cython.stochastic_swap'
+CYTHON_SOURCE_DIR = 'qiskit/transpiler/passes/mapping/cython/stochastic_swap'
+
+PACKAGE_DATA = {}
+
+INCLUDE_DIRS = []
+# Extra link args
+LINK_FLAGS = []
+# If on Win and not in MSYS2 (i.e. Visual studio compile)
+if (sys.platform == 'win32' and os.environ.get('MSYSTEM') is None):
+    COMPILER_FLAGS = ['/O2', '/std:c++11']
+# Everything else
+else:
+    COMPILER_FLAGS = ['-O2', '-funroll-loops', '-std=c++11']
+    if sys.platform == 'darwin':
+        # These are needed for compiling on OSX 10.14+
+        COMPILER_FLAGS.append('-mmacosx-version-min=10.9')
+        LINK_FLAGS.append('-mmacosx-version-min=10.9')
+
+# Remove -Wstrict-prototypes from cflags
+CFG_VARS = distutils.sysconfig.get_config_vars()
+if "CFLAGS" in CFG_VARS:
+    CFG_VARS["CFLAGS"] = CFG_VARS["CFLAGS"].replace("-Wstrict-prototypes", "")
+
+EXT_MODULES = []
+# Add Cython Extensions
+for ext in CYTHON_EXTS:
+    mod = Extension(CYTHON_MODULE+'.'+ext,
+                    sources=[CYTHON_SOURCE_DIR+'/'+ext+'.pyx'],
+                    include_dirs=INCLUDE_DIRS,
+                    extra_compile_args=COMPILER_FLAGS,
+                    extra_link_args=LINK_FLAGS,
+                    language='c++')
+    EXT_MODULES.append(mod)
+
 
 setup(
     name="qiskit-terra",
     version="0.8.0",
     description="Software for developing quantum computing programs",
-    long_description="""Terra provides the foundations for Qiskit. It allows the user to write 
+    long_description="""Terra provides the foundations for Qiskit. It allows the user to write
         quantum circuits easily, and takes care of the constraints of real hardware.""",
     url="https://github.com/Qiskit/qiskit-terra",
     author="Qiskit Development Team",
@@ -46,12 +88,16 @@ setup(
     ],
     keywords="qiskit sdk quantum",
     packages=find_packages(exclude=['test*']),
-    install_requires=requirements,
+    install_requires=REQUIREMENTS,
+    setup_requires=['Cython>=0.27.1'],
+    package_data=PACKAGE_DATA,
     include_package_data=True,
     python_requires=">=3.5",
     extra_requires={
         'visualization': ['matplotlib>=2.1', 'nxpd>=0.2', 'ipywidgets>=7.3.0',
                           'pydot'],
         'full-featured-simulators': ['qiskit-aer>=0.1']
-    }
+    },
+    ext_modules=cythonize(EXT_MODULES),
+    zip_safe=False
 )

test/python/transpiler/test_stochastic_swap.py

@@ -354,22 +354,22 @@ class TestStochasticSwap(QiskitTestCase):
         #  c_3: 0   ═══════════════════════════════════════════════════════╩═
         #
         expected = QuantumCircuit(qr, cr)
-        expected.x(qr[0])
-        expected.y(qr[1])
         expected.z(qr[2])
+        expected.y(qr[1])
+        expected.x(qr[0])
         expected.swap(qr[1], qr[2])
         expected.cx(qr[0], qr[2])
-        expected.measure(qr[0], cr[0])
         expected.swap(qr[2], qr[3])
-        expected.s(qr[3])
         expected.cx(qr[1], qr[2])
+        expected.s(qr[3])
         expected.t(qr[1])
         expected.h(qr[2])
-        expected.swap(qr[2], qr[3])
-        expected.cx(qr[2], qr[1])
-        expected.measure(qr[1], cr[2])
-        expected.measure(qr[2], cr[1])
-        expected.measure(qr[3], cr[3])
+        expected.measure(qr[0], cr[0])
+        expected.swap(qr[1], qr[2])
+        expected.cx(qr[3], qr[2])
+        expected.measure(qr[1], cr[3])
+        expected.measure(qr[3], cr[1])
+        expected.measure(qr[2], cr[2])
         expected_dag = circuit_to_dag(expected)
         #                      ┌───┐     ┌─┐
         # q_0: |0>─────────────┤ X ├──■──┤M├────────────────────────────────────────
@@ -514,20 +514,21 @@ class TestStochasticSwap(QiskitTestCase):
         # 0 - 1 - 3
         expected = QuantumCircuit(qr, ar, cr)
         expected.cx(qr[1], ar[0])
-        expected.h(ar[0])
-        expected.swap(qr[1], ar[1])
-        expected.cx(qr[0], qr[1])
+        expected.swap(qr[0], qr[1])
+        expected.cx(qr[1], ar[1])
         expected.h(ar[1])
-        expected.h(qr[1])
-        expected.cx(ar[0], qr[1])
-        expected.measure(qr[0], cr[0])
+        expected.h(ar[0])
+        expected.measure(qr[1], cr[0])
         expected.h(qr[0])
+        expected.swap(qr[1], ar[1])
+        expected.h(ar[1])
+        expected.cx(ar[0], qr[1])
         expected.measure(ar[0], cr[2])
         expected.swap(qr[1], ar[1])
-        expected.cx(qr[0], qr[1])
         expected.measure(ar[1], cr[3])
-        expected.measure(qr[1], cr[1])
-        expected.measure(qr[0], cr[0])
+        expected.cx(qr[1], qr[0])
+        expected.measure(qr[1], cr[0])
+        expected.measure(qr[0], cr[1])
         expected_dag = circuit_to_dag(expected)
 
         layout = Layout([(QuantumRegister(2, 'q'), 0),