Fixed GateMerge, KAK, gate def

- CZ gate def: required bits = 2.

- KAK: fixed numerical errors related to tolerance, added validation checks for simplified circuits.

- Gate merge: fixed endians, added validation.
Signed-off-by: Thien Nguyen <nguyentm@ornl.gov>

quantum/gate/ir/CommonGates.hpp

@@ -372,7 +372,7 @@ public:
   CZ(std::size_t srcqbit, std::size_t tgtqbit)
       : CZ(std::vector<std::size_t>{srcqbit, tgtqbit}) {}
 
-  const int nRequiredBits() const override { return 1; }
+  const int nRequiredBits() const override { return 2; }
 
   DEFINE_CLONE(CZ)
   DEFINE_VISITABLE()

quantum/plugins/circuits/kak/kak.cpp

@@ -197,9 +197,10 @@ bool isCanonicalized(double x, double y, double z)
 {
   // 0 ≤ abs(z) ≤ y ≤ x ≤ pi/4
   // if x = pi/4, z >= 0
-  if (std::abs(z) >= 0 && y >= std::abs(z) && x >= y && x <= M_PI_4)
+  const double TOL = 1e-9;
+  if (std::abs(z) >= 0 && y >= std::abs(z) && x >= y && x <= M_PI_4 + TOL)
   {
-    if (std::abs(x - M_PI_4) < 1e-9)
+    if (std::abs(x - M_PI_4) < TOL)
     {
       return (z >= 0);
     }
@@ -782,7 +783,72 @@ std::shared_ptr<CompositeInstruction> KAK::KakDecomposition::toGates(size_t in_b
       composite->addInstruction(gateRegistry->createInstruction("H", { bit1 }));
       composite->addInstruction(gateRegistry->createInstruction("Rx", { bit2 }, { -M_PI_2 }));
 
-      // TODO: add validation
+      const auto validateGateSequence = [&](const Eigen::Matrix4cd& in_target){
+        const auto H = []() {
+          GateMatrix result;
+          result << 1.0/std::sqrt(2), 1.0/std::sqrt(2), 1.0/std::sqrt(2), -1.0/std::sqrt(2);
+          return result;
+        };
+        const auto Rx = [](double angle) {
+          GateMatrix result;
+          result << std::cos(angle/2.0), -I*std::sin(angle/2.0), -I*std::sin(angle/2.0), std::cos(angle/2.0);
+          return result;
+        };
+        const auto Ry = [](double angle) {
+          GateMatrix result;
+          result << std::cos(angle/2), -std::sin(angle/2), std::sin(angle/2), std::cos(angle/2);
+          return result;
+        };
+        const auto Rz = [](double angle) {
+          GateMatrix result;
+          result << std::exp(-I*angle/2.0), 0, 0, std::exp(I*angle/2.0);
+          return result;
+        };
+        const auto CZ = []() {
+          Eigen::Matrix4cd cz;
+          cz << 1, 0, 0, 0, 
+                0, 1, 0, 0,
+                0, 0, 1, 0,
+                0, 0, 0, -1;
+          return cz;
+        };
+        
+        Eigen::Matrix2cd IdMat = Eigen::Matrix2cd::Identity();
+        Eigen::Matrix4cd totalU = Eigen::Matrix4cd::Identity();
+        totalU *= Eigen::kroneckerProduct(IdMat, Rx(-M_PI_2));
+        totalU *= Eigen::kroneckerProduct(H(), IdMat);
+        totalU *= CZ();
+        totalU *= Eigen::kroneckerProduct(H(), IdMat);
+        totalU *= Eigen::kroneckerProduct(IdMat, Rx(xAngle));
+        totalU *= Eigen::kroneckerProduct(Ry(yAngle), IdMat);
+        totalU *= Eigen::kroneckerProduct(H(), IdMat);
+        totalU *= CZ();
+        totalU *= Eigen::kroneckerProduct(H(), IdMat);      
+        totalU *= Eigen::kroneckerProduct(IdMat, Rx(M_PI_2));
+
+        // Find index of the largest element:
+        size_t colIdx = 0;
+        size_t rowIdx = 0;
+        double maxVal = std::abs(totalU(0,0));
+        for (size_t i = 0; i < totalU.rows(); ++i)
+        {
+          for (size_t j = 0; j < totalU.cols(); ++j)
+          {
+            if (std::abs(totalU(i,j)) > maxVal)
+            {
+              maxVal = std::abs(totalU(i,j));
+              colIdx = j;
+              rowIdx = i;
+            }
+          }
+        }
+
+        const std::complex<double> globalFactor = in_target(rowIdx, colIdx) / totalU(rowIdx, colIdx);
+        totalU = globalFactor * totalU;
+        return allClose(totalU, in_target);
+      };
+      
+      assert(validateGateSequence(interactionMatrixExp(x, y, z)));
       return composite;
     }
     // only XX is significant
@@ -796,7 +862,68 @@ std::shared_ptr<CompositeInstruction> KAK::KakDecomposition::toGates(size_t in_b
       composite->addInstruction(gateRegistry->createInstruction("CZ", { bit1, bit2 }));
       composite->addInstruction(gateRegistry->createInstruction("H", { bit1 }));
       
-      // TODO: add validation
+      const auto validateGateSequence = [&](const Eigen::Matrix4cd& in_target){
+        const auto H = []() {
+          GateMatrix result;
+          result << 1.0/std::sqrt(2), 1.0/std::sqrt(2), 1.0/std::sqrt(2), -1.0/std::sqrt(2);
+          return result;
+        };
+        const auto Rx = [](double angle) {
+          GateMatrix result;
+          result << std::cos(angle/2.0), -I*std::sin(angle/2.0), -I*std::sin(angle/2.0), std::cos(angle/2.0);
+          return result;
+        };
+        const auto Ry = [](double angle) {
+          GateMatrix result;
+          result << std::cos(angle/2), -std::sin(angle/2), std::sin(angle/2), std::cos(angle/2);
+          return result;
+        };
+        const auto Rz = [](double angle) {
+          GateMatrix result;
+          result << std::exp(-I*angle/2.0), 0, 0, std::exp(I*angle/2.0);
+          return result;
+        };
+        const auto CZ = []() {
+          Eigen::Matrix4cd cz;
+          cz << 1, 0, 0, 0, 
+                0, 1, 0, 0,
+                0, 0, 1, 0,
+                0, 0, 0, -1;
+          return cz;
+        };
+        
+        Eigen::Matrix2cd IdMat = Eigen::Matrix2cd::Identity();
+        Eigen::Matrix4cd totalU = Eigen::Matrix4cd::Identity();
+        
+        totalU *= Eigen::kroneckerProduct(H(), IdMat);
+        totalU *= CZ();
+        totalU *= Eigen::kroneckerProduct(IdMat, Rx(xAngle));
+        totalU *= CZ();
+        totalU *= Eigen::kroneckerProduct(H(), IdMat);      
+
+        // Find index of the largest element:
+        size_t colIdx = 0;
+        size_t rowIdx = 0;
+        double maxVal = std::abs(totalU(0,0));
+        for (size_t i = 0; i < totalU.rows(); ++i)
+        {
+          for (size_t j = 0; j < totalU.cols(); ++j)
+          {
+            if (std::abs(totalU(i,j)) > maxVal)
+            {
+              maxVal = std::abs(totalU(i,j));
+              colIdx = j;
+              rowIdx = i;
+            }
+          }
+        }
+
+        const std::complex<double> globalFactor = in_target(rowIdx, colIdx) / totalU(rowIdx, colIdx);
+        totalU = globalFactor * totalU;
+        return allClose(totalU, in_target);
+      };
+      
+      assert(validateGateSequence(interactionMatrixExp(x, y, z)));
       return composite; 
     }
   };

quantum/plugins/optimizers/gate_merge/GateMergeOptimizer.cpp

@@ -3,6 +3,53 @@
 #include "GateFusion.hpp"
 #include "xacc_service.hpp"
 
+namespace {
+bool compareMatIgnoreGlobalPhase(const Eigen::Matrix4cd& in_a, const Eigen::Matrix4cd& in_b)
+{
+    const auto flipKronOrder = [](const Eigen::Matrix4cd& in_mat){
+        Eigen::Matrix4cd result = Eigen::Matrix4cd::Zero();
+        const std::vector<size_t> order { 0, 2, 1, 3 };
+        for (size_t i = 0; i < in_mat.rows(); ++i)
+        {
+            for (size_t j = 0; j < in_mat.cols(); ++j)
+            {
+                result(order[i], order[j]) = in_mat(i, j);
+            }
+        }
+        return result;
+    };
+    
+    const auto bFixed = flipKronOrder(in_b);
+
+    // Find index of the largest element to normalize global phase.
+    size_t colIdx = 0;
+    size_t rowIdx = 0;
+    double maxVal = std::abs(in_a(0,0));
+    for (size_t i = 0; i < in_a.rows(); ++i)
+    {
+        for (size_t j = 0; j < in_a.cols(); ++j)
+        {
+            if (std::abs(in_a(i,j)) > maxVal)
+            {
+                maxVal = std::abs(in_a(i,j));
+                colIdx = j;
+                rowIdx = i;
+            }
+        }
+    }
+
+    const double TOL = 1e-6;
+    if (std::abs(std::abs(in_a(rowIdx, colIdx)) - std::abs(bFixed(rowIdx, colIdx))) > TOL)
+    {
+        return false;
+    }
+
+    const std::complex<double> globalFactor = in_a(rowIdx, colIdx) / bFixed(rowIdx, colIdx);
+    auto bFixedPhase = globalFactor * bFixed;
+    const auto diff = (bFixedPhase - in_a).norm();
+    return std::abs(diff) < TOL;
+}
+}
 namespace xacc {
 namespace quantum {
 void MergeSingleQubitGatesOptimizer::apply(std::shared_ptr<CompositeInstruction> program, const std::shared_ptr<Accelerator> accelerator, const HeterogeneousMap &options)
@@ -141,11 +188,11 @@ void MergeTwoQubitBlockOptimizer::apply(std::shared_ptr<CompositeInstruction> pr
                     assert(qubitPair.first  != qubitPair.second);
                     if (qubitPair.first < qubitPair.second)
                     {
-                        return (bit == qubitPair.first) ? 0 : 1;
+                        return (bit == qubitPair.first) ? 1 : 0;
                     }
                     else
                     {
-                        return (bit == qubitPair.first) ? 1 : 0;
+                        return (bit == qubitPair.first) ? 0 : 1;
                     }
                 }; 
 
@@ -194,11 +241,17 @@ void MergeTwoQubitBlockOptimizer::apply(std::shared_ptr<CompositeInstruction> pr
                 std::make_pair("unitary", uMat)
             });
             assert(expandOk);
-
+            
+            const auto calcUopt = [](const std::shared_ptr<CompositeInstruction> composite) {
+                auto fuser = xacc::getService<xacc::quantum::GateFuser>("default");
+                fuser->initialize(composite);
+                return fuser->calcFusedGate(2);
+            };
+            
             // Optimized decomposed sequence:
             const auto nbInstructionsAfter = kak->nInstructions();
             // A simplified sequence was found.
-            if (nbInstructionsAfter < sequence.size())
+            if (nbInstructionsAfter < sequence.size() && compareMatIgnoreGlobalPhase(uMat, calcUopt(kak)))
             {
                 // Disable to remove:
                 const auto programLengthBefore = program->nInstructions();

quantum/plugins/optimizers/gate_merge/tests/GateMergingTester.cpp

@@ -72,32 +72,40 @@ TEST(GateMergingTester, checkTwoQubitSimple)
     auto c = xacc::getService<xacc::Compiler>("xasm");
     auto f = c->compile(R"(__qpu__ void test3(qbit q) {
         H(q[0]);
-        Z(q[2]);
-        CNOT(q[2], q[1]);
-        H(q[2]);
-        T(q[1]);
-        X(q[0]);
-        CNOT(q[1], q[2]);
-        H(q[2]);
+        H(q[1]);
+        Rz(q[0], 1.234);
         Y(q[1]);
-        CNOT(q[3], q[4]);
-        X(q[2]);
-        X(q[1]);
-        CNOT(q[2], q[1]);
-        H(q[2]);
-        T(q[1]);
+        T(q[0]);
+        Z(q[1]);
+        H(q[0]);
+        H(q[1]);
+        S(q[0]);
+        H(q[1]);
         X(q[0]);
-        CNOT(q[1], q[2]);
-        H(q[2]);
-        Y(q[1]);
-        CNOT(q[3], q[4]);
-        X(q[2]);
-        X(q[1]);
+        Z(q[1]);
         CNOT(q[1], q[0]);
         H(q[0]);
         H(q[1]);
+        CPhase(q[0], q[1], 1.123);
+        H(q[0]);
+        T(q[1]);
+        X(q[0]);
+        CNOT(q[0], q[1]);
+        T(q[0]);
+        //H(q[2]);
+        Rx(q[0], 1.234);
+        Ry(q[1], -2.456);
+        H(q[1]);
+        Rz(q[0], 3.124);
+        H(q[0]);
+        H(q[1]);
+        Y(q[0]);
+        H(q[0]);
+        H(q[1]);
+        T(q[0]);
+        X(q[1]);
     })")->getComposites()[0];
-
+   
     auto opt = xacc::getService<xacc::IRTransformation>("two-qubit-block-merging");
     const auto nbInstBefore = f->nInstructions();
     auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
@@ -107,19 +115,15 @@ TEST(GateMergingTester, checkTwoQubitSimple)
         circuitCopy->addInstruction(f->getInstruction(i)->clone());
     }
     opt->apply(f, nullptr);
-    const auto nbInstAfter = f->nInstructions();
-    std::cout << "Before: " << nbInstBefore << "; After: " <<  nbInstAfter << "\n";
     std::cout << "HOWDY:\n" << f->toString() << "\n";
-    EXPECT_TRUE(nbInstAfter < nbInstBefore);
-    EXPECT_TRUE(circuitCopy->nInstructions() == nbInstBefore);
-    std::cout << "HOWDY:\n" << circuitCopy->toString() << "\n";
-
+    const auto nbInstAfter = f->nInstructions();
+    EXPECT_TRUE(nbInstAfter < nbInstBefore); 
     // Validate using gate fusion:
     auto fuser = xacc::getService<xacc::quantum::GateFuser>("default");
     fuser->initialize(circuitCopy);
-    const Eigen::MatrixXcd uMatOriginal = fuser->calcFusedGate(5);
+    const Eigen::MatrixXcd uMatOriginal = fuser->calcFusedGate(2);
     fuser->initialize(f);
-    const Eigen::MatrixXcd uMatAfter = fuser->calcFusedGate(5);
+    const Eigen::MatrixXcd uMatAfter = fuser->calcFusedGate(2);
     // Compensate any global phase differences.
     // Find index of the largest element:
     size_t colIdx = 0;