Connectivity constraints for clifford synthesis

include/cliffordsynthesis/CliffordSynthesizer.hpp

@@ -12,8 +12,10 @@
 #include "cliffordsynthesis/TargetMetric.hpp"
 #include "cliffordsynthesis/encoding/SATEncoder.hpp"
 #include "ir/QuantumComputation.hpp"
+#include "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

@@ -16,6 +16,8 @@
 #include <ostream>
 #include <sstream>
 #include <string>
+#include <utility>
+#include <vector>
 
 namespace cs {
 class Results {
@@ -53,6 +55,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 : pvector) {
+      for (const bool val : row) {
+        oss << (val ? '1' : '0');
+      }
+      oss << '\n';
+    }
+    return oss.str();
+  }
+  [[nodiscard]] std::vector<std::vector<bool>> getMappingVector() const {
+    return pvector;
+  }
+
   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; }
@@ -70,6 +86,17 @@ class Results {
     ss << tableau;
     resultTableau = ss.str();
   }
+  void setMapping(std::vector<std::vector<bool>> p) {
+    std::ostringstream oss;
+    for (const auto& row : pvector) {
+      for (const bool val : row) {
+        oss << (val ? '1' : '0');
+      }
+      oss << '\n';
+    }
+    pvals = oss.str();
+    pvector = std::move(p);
+  }
 
   [[nodiscard]] bool sat() const {
     return getSolverResult() == logicbase::Result::SAT;
@@ -103,6 +130,8 @@ class Results {
   double runtime = 0.0;
   std::size_t solverCalls = 0U;
 
+  std::string pvals;
+  std::vector<std::vector<bool>> pvector;
   std::string resultTableau;
   std::string resultCircuit;
 };

include/cliffordsynthesis/Tableau.hpp

@@ -163,5 +163,11 @@ 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 reverseMapping(const std::vector<std::vector<bool>>* p);
+  Tableau reverseMappingOnRows(const std::vector<std::vector<bool>>* p,
+                               size_t nq);
+  void gaussianEliminationGF2();
+  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 "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 "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/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;
@@ -491,6 +493,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;

src/cliffordsynthesis/Tableau.cpp

@@ -492,4 +492,107 @@ bool Tableau::isIdentityTableau() const {
   }
   return true;
 }
+
+Tableau Tableau::applyMapping(const std::vector<std::vector<bool>>* p) {
+  Tableau mappedTableau = Tableau(nQubits, hasDestabilizers());
+  for (size_t i = 0; i < mappedTableau.getTableauSize(); i++) {
+    for (unsigned char& j : mappedTableau.tableau[i]) {
+      j = 0;
+    }
+  }
+  for (size_t i = 0; i < p->size(); i++) {
+    for (size_t j = 0; j < (*p)[i].size(); j++) {
+      // apply mapping from column i to j if p is set
+      if ((*p)[i][j]) {
+        // in every row swap x entry and z entry
+        for (size_t n = 0; n < mappedTableau.getTableauSize(); n++) {
+          mappedTableau.tableau[n][j] = tableau[n][i];
+          mappedTableau.tableau[n][j + mappedTableau.nQubits] =
+              tableau[n][i + mappedTableau.nQubits];
+        }
+      }
+    }
+  }
+  // copy r column without changes
+  for (size_t i = 0; i < tableau.size(); i++) {
+    mappedTableau.tableau[i][2 * nQubits] = tableau[i][2 * nQubits];
+  }
+  return mappedTableau;
+}
+
+// number of Qubits is passed because nQubits is not set in result Tableau of
+// synthesis
+Tableau Tableau::reverseMappingOnRows(const std::vector<std::vector<bool>>* p,
+                                      const size_t nq) {
+  auto mappedTableau = Tableau(nq, true);
+  mappedTableau.tableau = tableau;
+  for (size_t i = 0; i < p->size(); i++) {
+    for (size_t j = 0; j < (*p)[i].size(); j++) {
+      // apply mapping from row i to j if p is set
+      if ((*p)[i][j]) {
+        mappedTableau.tableau[i] = tableau[j];
+        mappedTableau.tableau[i + nq] = tableau[j + nq];
+      }
+    }
+  }
+  return mappedTableau;
+}
+
+// in place Gauss Elimination
+void Tableau::gaussianEliminationGF2() {
+  const size_t rows = tableau.size();
+  const size_t cols = tableau[0].size();
+  if (rows == 1) {
+    return;
+  }
+
+  size_t pivot = 0;
+
+  for (size_t col = 0; col < cols; ++col) {
+    // find the pivot row for the current column
+    size_t maxRow = pivot;
+    if (maxRow >= rows) {
+      break;
+    }
+    for (size_t row = pivot; row < rows; ++row) {
+      if (tableau[row][col] == 1) {
+        maxRow = row;
+        break;
+      }
+    }
+    // if no pivot is found, continue to the next column
+    if (tableau[maxRow][col] == 0) {
+      continue;
+    }
+    // swap the pivot row with the current row
+    std::swap(tableau[pivot], tableau[maxRow]);
+
+    // eliminate all other 1s in the current column
+    for (size_t row = 0; row < rows; ++row) {
+      if (row != pivot && tableau[row][col] == 1) {
+        for (size_t k = 0; k < cols; ++k) {
+          if (tableau[row][k] == tableau[pivot][k]) {
+            tableau[row][k] = 0;
+          } else {
+            tableau[row][k] = 1;
+          }
+        }
+      }
+    }
+    pivot++;
+  }
+}
+// equivalence check for two Tableaus without destablilizers
+bool Tableau::equivalentUpToStabilizer(const Tableau* t) const {
+  assert(getTableauSize() == t->getTableauSize());
+
+  Tableau tEliminated = Tableau(nQubits, hasDestabilizers());
+  tEliminated.tableau = tableau;
+  tEliminated.gaussianEliminationGF2();
+
+  Tableau t2Eliminated = Tableau(nQubits, hasDestabilizers());
+  t2Eliminated.tableau = t->getTableau();
+  t2Eliminated.gaussianEliminationGF2();
+  return tEliminated == t2Eliminated;
+}
 } // namespace cs