Connectivity constraints for clifford synthesis

docs/source/Synthesis.ipynb

@@ -286,6 +286,55 @@
     "qc_synth.draw(output=\"mpl\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Connectivity Constraints on qubits"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "Limited qubit connectivity is implemented to handle the physical constraints of quantum hardware, where certain qubits can only interact with specific other qubits.\n",
+    "For this a custom coupling map can be passed to the synthesizer, which ensures that the generated circuit from the synthesis respects these connectivity constraints.\n",
+    "The mapping used in the synthesis is also stored in the result returned by the synthesis\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit import QuantumCircuit\n",
+    "\n",
+    "from mqt import qmap\n",
+    "\n",
+    "circ = QuantumCircuit(3)\n",
+    "circ.h(1)\n",
+    "circ.h(0)\n",
+    "circ.cx(0, 1)\n",
+    "circ.cx(0, 2)\n",
+    "circ.draw(output=\"mpl\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "circ_opt, results = qmap.optimize_clifford(circ, coupling_map={(1, 0), (0, 1), (1, 2), (2, 1)})\n",
+    "\n",
+    "print(results.mapping)\n",
+    "print(results.mapping_string)\n",
+    "circ_opt.draw(output=\"mpl\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

include/cliffordsynthesis/CliffordSynthesizer.hpp

@@ -12,8 +12,10 @@
 #include "cliffordsynthesis/TargetMetric.hpp"
 #include "cliffordsynthesis/encoding/SATEncoder.hpp"
 #include "ir/QuantumComputation.hpp"
+#include "sc/utils.hpp"
 
 #include <cstddef>
+#include <cstdint>
 #include <limits>
 #include <memory>
 #include <plog/Log.h>
@@ -28,23 +30,49 @@ class CliffordSynthesizer {
 public:
   CliffordSynthesizer() = default;
   CliffordSynthesizer(Tableau initial, Tableau target)
-      : initialTableau(std::move(initial)), targetTableau(std::move(target)) {}
+      : initialTableau(std::move(initial)),
+        couplingMap(getFullyConnectedMap(
+            static_cast<uint16_t>(target.getQubitCount()))),
+        targetTableau(std::move(target)) {}
+  CliffordSynthesizer(Tableau initial, Tableau target, CouplingMap qm)
+      : initialTableau(std::move(initial)), couplingMap(std::move(qm)),
+        targetTableau(std::move(target)) {}
   explicit CliffordSynthesizer(Tableau target)
       : initialTableau(target.getQubitCount(), target.hasDestabilizers()),
+        couplingMap(getFullyConnectedMap(
+            static_cast<uint16_t>(target.getQubitCount()))),
         targetTableau(std::move(target)) {}
   CliffordSynthesizer(Tableau initial, qc::QuantumComputation& qc)
       : initialTableau(std::move(initial)),
+        couplingMap(
+            getFullyConnectedMap(static_cast<uint16_t>(qc.getNqubits()))),
         targetTableau(qc, 0, std::numeric_limits<std::size_t>::max(),
                       initialTableau.hasDestabilizers()),
-        initialCircuit(std::make_shared<qc::QuantumComputation>(qc)),
-        results(qc, targetTableau) {}
+        initialCircuit(std::make_shared<qc::QuantumComputation>(qc)) {}
   explicit CliffordSynthesizer(qc::QuantumComputation& qc,
                                const bool useDestabilizers = false)
       : initialTableau(qc.getNqubits(), useDestabilizers),
+        couplingMap(
+            getFullyConnectedMap(static_cast<uint16_t>(qc.getNqubits()))),
+        targetTableau(qc, 0, std::numeric_limits<std::size_t>::max(),
+                      useDestabilizers),
+        initialCircuit(std::make_shared<qc::QuantumComputation>(qc)) {}
+  explicit CliffordSynthesizer(Tableau target, CouplingMap qm)
+      : initialTableau(target.getQubitCount(), target.hasDestabilizers()),
+        couplingMap(std::move(qm)), targetTableau(std::move(target)) {}
+  CliffordSynthesizer(Tableau initial, qc::QuantumComputation& qc,
+                      CouplingMap qm)
+      : initialTableau(std::move(initial)), couplingMap(std::move(qm)),
+        targetTableau(qc, 0, std::numeric_limits<std::size_t>::max(),
+                      initialTableau.hasDestabilizers()),
+        initialCircuit(std::make_shared<qc::QuantumComputation>(qc)) {}
+  explicit CliffordSynthesizer(qc::QuantumComputation& qc, CouplingMap qm,
+                               const bool useDestabilizers = false)
+      : initialTableau(qc.getNqubits(), useDestabilizers),
+        couplingMap(std::move(qm)),
         targetTableau(qc, 0, std::numeric_limits<std::size_t>::max(),
                       useDestabilizers),
-        initialCircuit(std::make_shared<qc::QuantumComputation>(qc)),
-        results(qc, targetTableau) {}
+        initialCircuit(std::make_shared<qc::QuantumComputation>(qc)) {}
 
   virtual ~CliffordSynthesizer() = default;
 
@@ -72,6 +100,7 @@ class CliffordSynthesizer {
 
 protected:
   Tableau initialTableau;
+  CouplingMap couplingMap;
   Tableau targetTableau;
   std::shared_ptr<qc::QuantumComputation> initialCircuit;
 

include/cliffordsynthesis/Results.hpp

@@ -13,7 +13,10 @@
 #include <limits>
 #include <nlohmann/json_fwd.hpp>
 #include <ostream>
+#include <sstream>
 #include <string>
+#include <utility>
+#include <vector>
 
 namespace cs {
 class Results {
@@ -40,6 +43,20 @@ class Results {
   [[nodiscard]] std::string getResultCircuit() const { return resultCircuit; }
   [[nodiscard]] std::string getResultTableau() const { return resultTableau; }
 
+  [[nodiscard]] std::string getMapping() const {
+    std::ostringstream oss;
+    for (const auto& row : permutationVector) {
+      for (const bool val : row) {
+        oss << (val ? '1' : '0');
+      }
+      oss << '\n';
+    }
+    return oss.str();
+  }
+  [[nodiscard]] std::vector<std::vector<bool>> getMappingVector() const {
+    return permutationVector;
+  }
+
   void setSingleQubitGates(const std::size_t g) { singleQubitGates = g; }
   void setTwoQubitGates(const std::size_t g) { twoQubitGates = g; }
   void setDepth(const std::size_t d) { depth = d; }
@@ -49,6 +66,17 @@ class Results {
 
   void setResultCircuit(qc::QuantumComputation& qc);
   void setResultTableau(const Tableau& tableau);
+  void setMapping(std::vector<std::vector<bool>> p) {
+    std::ostringstream oss;
+    for (const auto& row : permutationVector) {
+      for (const bool val : row) {
+        oss << (val ? '1' : '0');
+      }
+      oss << '\n';
+    }
+    permutationString = oss.str();
+    permutationVector = std::move(p);
+  }
 
   [[nodiscard]] bool sat() const {
     return getSolverResult() == logicbase::Result::SAT;
@@ -69,6 +97,8 @@ class Results {
   double runtime = 0.0;
   std::size_t solverCalls = 0U;
 
+  std::string permutationString;
+  std::vector<std::vector<bool>> permutationVector;
   std::string resultTableau;
   std::string resultCircuit;
 };

include/cliffordsynthesis/Tableau.hpp

@@ -163,5 +163,10 @@ class Tableau {
   [[nodiscard]] std::uint64_t getBVFrom(const std::size_t column) const {
     return getBVFrom<64>(column).to_ullong();
   }
+  Tableau applyMapping(const std::vector<std::vector<bool>>& p);
+  Tableau reverseMappingOnRows(const std::vector<std::vector<bool>>& p,
+                               size_t nq);
+  void rref();
+  bool equivalentUpToStabilizer(const Tableau* t) const;
 };
 } // namespace cs

include/cliffordsynthesis/encoding/GateEncoder.hpp

@@ -12,6 +12,7 @@
 #include "logicblocks/Logic.hpp"
 #include "logicblocks/LogicBlock.hpp"
 #include "logicblocks/LogicTerm.hpp"
+#include "sc/utils.hpp"
 
 #include <array>
 #include <cstddef>
@@ -27,8 +28,9 @@ class GateEncoder {
   GateEncoder(const std::size_t nQubits, const std::size_t tableauSize,
               const std::size_t timestepLimit,
               TableauEncoder::Variables* tableauVars,
-              std::shared_ptr<logicbase::LogicBlock> logicBlock)
-      : N(nQubits), S(tableauSize), T(timestepLimit), tvars(tableauVars),
+              std::shared_ptr<logicbase::LogicBlock> logicBlock, CouplingMap cm)
+      : N(nQubits), S(tableauSize), T(timestepLimit),
+        couplingMap(std::move(cm)), tvars(tableauVars),
         lb(std::move(logicBlock)) {}
   virtual ~GateEncoder() = default;
 
@@ -115,6 +117,9 @@ class GateEncoder {
   // timestep limit T
   std::size_t T{}; // NOLINT (readability-identifier-naming)
 
+  // coupling Map
+  CouplingMap couplingMap;
+
   // the gate variables
   Variables vars{};
 

include/cliffordsynthesis/encoding/SATEncoder.hpp

@@ -14,6 +14,7 @@
 #include "cliffordsynthesis/encoding/TableauEncoder.hpp"
 #include "logicblocks/Logic.hpp"
 #include "logicblocks/LogicBlock.hpp"
+#include "sc/utils.hpp"
 
 #include <cstddef>
 #include <memory>
@@ -33,6 +34,9 @@ class SATEncoder {
     // the number of qubits to encode
     std::size_t nQubits{};
 
+    // coupling Map of Qubits
+    CouplingMap couplingMap;
+
     // the number of timesteps to encode
     std::size_t timestepLimit{};
 

include/cliffordsynthesis/encoding/TableauEncoder.hpp

@@ -33,6 +33,8 @@ class TableauEncoder {
     logicbase::LogicMatrix z;
     // variables for the phase parts of the tableaus
     logicbase::LogicVector r;
+    // variables for mapping of qubits
+    logicbase::LogicMatrix p;
 
     // update rules for single-qubit gates
     [[nodiscard]] logicbase::LogicTerm
@@ -58,14 +60,20 @@ class TableauEncoder {
   void createTableauVariables();
 
   // fixing the tableau
-  void assertTableau(const Tableau& tableau, std::size_t t);
+  void assertTableau(const Tableau& tableau, std::size_t t) const;
 
   // extracting the tableau
   void extractTableauFromModel(Results& results, std::size_t t,
                                logicbase::Model& model) const;
 
   [[nodiscard]] auto* getVariables() { return &vars; }
 
+  // get mapping variables and store them in results
+  void extractMappingFromModel(Results& results, logicbase::Model& model) const;
+
+  // assert constraints for mapping variables
+  void assertMappingConstraints() const;
+
 protected:
   // number of qubits N
   std::size_t N{}; // NOLINT (readability-identifier-naming)

src/cliffordsynthesis/CMakeLists.txt

@@ -10,11 +10,12 @@ if(NOT TARGET ${MQT_QMAP_CLIFFORD_SYNTHESIS_TARGET_NAME})
   file(GLOB_RECURSE CLIFFORD_SYNTHESIS_HEADERS
        ${MQT_QMAP_INCLUDE_BUILD_DIR}/cliffordsynthesis/*.hpp)
   file(GLOB_RECURSE CLIFFORD_SYNTHESIS_SOURCES **.cpp)
+  list(APPEND CLIFFORD_SYNTHESIS_SOURCES "${CMAKE_SOURCE_DIR}/src/sc/utils.cpp")
 
   # add CliffordSynthesis Package library
   add_library(${MQT_QMAP_CLIFFORD_SYNTHESIS_TARGET_NAME} ${CLIFFORD_SYNTHESIS_HEADERS}
                                                          ${CLIFFORD_SYNTHESIS_SOURCES})
-
+  message(${CLIFFORD_SYNTHESIS_SOURCES})
   # set include directories
   target_include_directories(${MQT_QMAP_CLIFFORD_SYNTHESIS_TARGET_NAME}
                              PUBLIC $<BUILD_INTERFACE:${MQT_QMAP_INCLUDE_BUILD_DIR}>)

src/cliffordsynthesis/CliffordSynthesizer.cpp

@@ -6,6 +6,7 @@
 #include "cliffordsynthesis/CliffordSynthesizer.hpp"
 
 #include "cliffordsynthesis/Configuration.hpp"
+#include "cliffordsynthesis/Results.hpp"
 #include "cliffordsynthesis/Tableau.hpp"
 #include "cliffordsynthesis/TargetMetric.hpp"
 #include "cliffordsynthesis/encoding/SATEncoder.hpp"
@@ -52,6 +53,7 @@ void CliffordSynthesizer::synthesize(const Configuration& config) {
   encoderConfig.initialTableau = &initialTableau;
   encoderConfig.targetTableau = &targetTableau;
   encoderConfig.nQubits = initialTableau.getQubitCount();
+  encoderConfig.couplingMap = couplingMap;
   encoderConfig.timestepLimit = configuration.initialTimestepLimit;
   encoderConfig.targetMetric = configuration.target;
   encoderConfig.useMaxSAT = configuration.useMaxSAT;
@@ -85,6 +87,7 @@ void CliffordSynthesizer::synthesize(const Configuration& config) {
     const auto [lowerBin, upperBin] = determineUpperBound(encoderConfig);
     lower = lowerBin;
     upper = upperBin;
+    PLOG_INFO << "Upper bound " << upperBin;
 
     // if the upper bound is 0, the solution does not require any gates and the
     // synthesis is done.
@@ -95,7 +98,10 @@ void CliffordSynthesizer::synthesize(const Configuration& config) {
   }
   // Otherwise, the determined upper bound is used as an initial timestep
   // limit.
+
   encoderConfig.timestepLimit = upper;
+  PLOG_INFO << "Upper bound " << encoderConfig.timestepLimit;
+
   // Once a valid upper bound is found, the SAT problem is solved again with
   // the objective function encoded.
   switch (config.target) {
@@ -491,6 +497,7 @@ std::vector<std::size_t> getLayers(const qc::QuantumComputation& qc) {
 void CliffordSynthesizer::depthHeuristicSynthesis() {
   PLOG_INFO << "Optimizing Circuit with Heuristic";
   if (initialCircuit->getDepth() == 0) {
+    results.setResultCircuit(*initialCircuit);
     return;
   }
   auto optimalConfig = configuration;