Implemented two-qubit block merging

Signed-off-by: Thien Nguyen <nguyentm@ornl.gov>

quantum/plugins/optimizers/gate_merge/GateMergeOptimizer.cpp

@@ -30,7 +30,7 @@ void MergeSingleQubitGatesOptimizer::apply(std::shared_ptr<CompositeInstruction>
             const bool expandOk = zyz->expand({ 
                 std::make_pair("unitary", uMat)
             });
-
+            assert(expandOk);
             // Optimized decomposed sequence:
             const auto nbInstructionsAfter = zyz->nInstructions();
             // A simplified sequence was found.
@@ -62,7 +62,6 @@ void MergeSingleQubitGatesOptimizer::apply(std::shared_ptr<CompositeInstruction>
 std::vector<size_t> MergeSingleQubitGatesOptimizer::findSingleQubitGateSequence(const std::shared_ptr<CompositeInstruction> in_program, size_t in_startIdx, size_t in_lengthLimit) const
 {
     const auto nbInstructions = in_program->nInstructions();
-    std::unordered_map<size_t, std::vector<size_t>> qubitToSequence;
     assert(in_startIdx < nbInstructions);
     auto firstInst = in_program->getInstruction(in_startIdx);
     // Not a single-qubit gate.
@@ -104,5 +103,200 @@ std::vector<size_t> MergeSingleQubitGatesOptimizer::findSingleQubitGateSequence(
     // Reach the end of the circuit:  
     return returnSeq(gateSequence);    
 }
+
+void MergeTwoQubitBlockOptimizer::apply(std::shared_ptr<CompositeInstruction> program, const std::shared_ptr<Accelerator> accelerator, const HeterogeneousMap& options)
+{
+    auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
+
+    // No need to optimize block with less than 6 gates
+    // since the KAK decomposition will result in at least 5 gate.
+    const size_t MIN_SIZE = 6;
+    if (program->nInstructions() < MIN_SIZE)
+    {
+        return;
+    }
+    for (size_t instIdx = 0; instIdx < program->nInstructions(); ++instIdx)
+    {
+        std::pair<size_t, size_t> qubitPair = std::make_pair(0, 0);
+        const auto sequence = findGateSequence(program, instIdx, MIN_SIZE, qubitPair);
+        if (!sequence.empty()) 
+        {
+            auto tmpKernel = gateRegistry->createComposite("__TMP__");
+            
+            const auto mapBits = [&qubitPair](const std::vector<size_t>& in_bits){
+                const auto translate = [&qubitPair](size_t bit) {
+                    assert(bit == qubitPair.first || bit == qubitPair.second);
+                    assert(qubitPair.first  != qubitPair.second);
+                    if (qubitPair.first < qubitPair.second)
+                    {
+                        return (bit == qubitPair.first) ? 0 : 1;
+                    }
+                    else
+                    {
+                        return (bit == qubitPair.first) ? 1 : 0;
+                    }
+                }; 
+
+                std::vector<size_t> newBits;
+                for (const auto& bit: in_bits)         
+                {
+                    newBits.emplace_back(translate(bit));
+                }   
+                return newBits;
+            };
+
+            // Map { 0, 1 } bits back to original bits
+            const auto remapBits = [&qubitPair](const std::vector<size_t>& in_bits){
+                const auto translate = [&qubitPair](size_t bit) {
+                    assert(bit == 0 || bit == 1);
+                    assert(qubitPair.first  != qubitPair.second);
+                    if (qubitPair.first < qubitPair.second)
+                    {
+                        return (bit == 0) ? qubitPair.first : qubitPair.second;
+                    }
+                    else
+                    {
+                        return (bit == 0) ? qubitPair.second : qubitPair.first;
+                    }
+                }; 
+
+                std::vector<size_t> newBits;
+                for (const auto& bit: in_bits)         
+                {
+                    newBits.emplace_back(translate(bit));
+                }   
+                return newBits;
+            };
+
+            for (const auto& instIdx: sequence)
+            {
+                auto instrPtr = program->getInstruction(instIdx)->clone();
+                instrPtr->setBits(mapBits(instrPtr->bits()));
+                tmpKernel->addInstruction(instrPtr);
+            }
+            auto fuser = xacc::getService<xacc::quantum::GateFuser>("default");
+            fuser->initialize(tmpKernel);
+            const Eigen::Matrix4cd uMat = fuser->calcFusedGate(2);
+            auto kak = std::dynamic_pointer_cast<quantum::Circuit>(xacc::getService<Instruction>("kak"));
+            const bool expandOk = kak->expand({ 
+                std::make_pair("unitary", uMat)
+            });
+            assert(expandOk);
+
+            // Optimized decomposed sequence:
+            const auto nbInstructionsAfter = kak->nInstructions();
+            // A simplified sequence was found.
+            if (nbInstructionsAfter < sequence.size())
+            {
+                // Disable to remove:
+                const auto programLengthBefore = program->nInstructions();
+                for (const auto& instIdx: sequence)
+                {
+                    auto instrPtr = program->getInstruction(instIdx);
+                    instrPtr->disable();
+                }
+                program->removeDisabled();
+
+                const auto locationToInsert = sequence[0];
+                for (auto& newInst: kak->getInstructions())
+                {
+                    newInst->setBits(remapBits(newInst->bits()));
+                    program->insertInstruction(locationToInsert, newInst->clone());
+                }
+                const auto programLengthAfter = program->nInstructions();
+                assert(programLengthAfter < programLengthBefore);
+            }
+        }
+    }
+}
+
+std::vector<size_t> MergeTwoQubitBlockOptimizer::findGateSequence(const std::shared_ptr<CompositeInstruction> in_program, size_t in_startIdx, size_t in_lengthLimit, std::pair<size_t, size_t>& out_qubitPair) const
+{
+    const auto nbInstructions = in_program->nInstructions();
+    assert(in_startIdx < nbInstructions);
+    auto firstInst = in_program->getInstruction(in_startIdx);
+    if (firstInst->name() == "Measure")
+    {
+        return {};
+    }
+
+    const auto qubitPairIfAny = [&]() -> std::optional<std::pair<size_t, size_t>> {
+        if (firstInst->bits().size() == 2)
+        {
+            return std::make_pair(firstInst->bits()[0], firstInst->bits()[1]);
+        }
+
+        // Single qubit gate:
+        // Scan forward to find the *first* two-qubit gate which involves this qubit.
+        assert(firstInst->bits().size() == 1);
+        const auto firstBitIdx =  firstInst->bits()[0];
+        for (size_t instIdx = in_startIdx + 1; instIdx < nbInstructions; ++instIdx)
+        {
+            auto instPtr = in_program->getInstruction(instIdx);
+            if (instPtr->bits().size() == 2 && xacc::container::contains(instPtr->bits(), firstBitIdx))
+            {
+                return std::make_pair(instPtr->bits()[0], instPtr->bits()[1]);
+            }
+        }
+        // Cannot find the boundary: i.e. no two-qubit gates
+        return std::optional<std::pair<size_t, size_t>>();
+    }();
+
+    if (!qubitPairIfAny.has_value())
+    {
+        return {};
+    }
+    const auto qubitPair = qubitPairIfAny.value();
+    std::vector<size_t> gateSequence;
+    gateSequence.emplace_back(in_startIdx);
+    out_qubitPair = qubitPair;
+    const auto returnSeq = [&](const std::vector<size_t>& in_seq) -> std::vector<size_t> {
+        return (in_seq.size() >= in_lengthLimit) ? in_seq : std::vector<size_t>{};
+    };
+
+    for (size_t instIdx = in_startIdx + 1; instIdx < nbInstructions; ++instIdx)
+    {
+        auto instPtr = in_program->getInstruction(instIdx);
+        if (instPtr->bits().size() == 1)
+        {
+            if (instPtr->bits()[0] == qubitPair.first || instPtr->bits()[0] == qubitPair.second)
+            {
+                if (instPtr->name() == "Measure")
+                {
+                    return returnSeq(gateSequence);
+                }
+                else
+                {
+                    gateSequence.emplace_back(instIdx);
+                }
+            }
+        }
+        else if (instPtr->bits().size() == 2)
+        {
+            const auto areQubitOperandsMatched = [&qubitPair](const std::vector<size_t>& in_bits) {
+                const bool match1 =  (in_bits[0] == qubitPair.first) && (in_bits[1] == qubitPair.second); 
+                const bool match2 =  (in_bits[1] == qubitPair.first) && (in_bits[0] == qubitPair.second); 
+                return match1 || match2;
+            };
+
+            const auto hasOnlyOneBitMatched = [&qubitPair](const std::vector<size_t>& in_bits) {
+                const bool oneMatched = xacc::container::contains(in_bits, qubitPair.first) || xacc::container::contains(in_bits, qubitPair.second);
+                const bool oneNotMatched = !xacc::container::contains(in_bits, qubitPair.first) || !xacc::container::contains(in_bits, qubitPair.second);
+                return oneMatched && oneNotMatched;
+            };
+
+            if (areQubitOperandsMatched(instPtr->bits()))
+            {
+                gateSequence.emplace_back(instIdx);
+            }
+            else if (hasOnlyOneBitMatched(instPtr->bits()))
+            {
+                return returnSeq(gateSequence);
+            }
+        }
+    }
+    // End of circuit:
+    return returnSeq(gateSequence);
+}
 }
 }
\ No newline at end of file

quantum/plugins/optimizers/gate_merge/GateMergeOptimizer.hpp

@@ -29,10 +29,14 @@ class MergeTwoQubitBlockOptimizer : public IRTransformation
 public:
     virtual void apply(std::shared_ptr<CompositeInstruction> program,
                      const std::shared_ptr<Accelerator> accelerator,
-                     const HeterogeneousMap &options = {}) override { /* TODO */ };
+                     const HeterogeneousMap &options = {}) override;
     virtual const IRTransformationType type() const override { return IRTransformationType::Optimization; }
     const std::string name() const override { return "two-qubit-block-merging"; }
     const std::string description() const override { return ""; }
+
+private:
+    // Finds the sequence of gates (from the start index) which forms a two-qubit block (with no connections to outside the block)
+    std::vector<size_t> findGateSequence(const std::shared_ptr<CompositeInstruction> in_program, size_t in_startIdx, size_t in_lengthLimit, std::pair<size_t, size_t>& out_qubitPair) const;    
 };
 }
 }
\ No newline at end of file

quantum/plugins/optimizers/gate_merge/tests/CMakeLists.txt

 add_xacc_test(GateMerging)
-target_link_libraries(GateMergingTester xacc-quantum-gate)
\ No newline at end of file
+target_link_libraries(GateMergingTester xacc-quantum-gate xacc-circuit-optimizers)
\ No newline at end of file

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

@@ -3,7 +3,7 @@
 #include "xacc.hpp"
 #include "xacc_service.hpp"
 #include "IRTransformation.hpp"
-
+#include "GateFusion.hpp"
 
 TEST(GateMergingTester, checkSingleQubitSimple) 
 {
@@ -67,6 +67,88 @@ TEST(GateMergingTester, checkMixing)
     EXPECT_EQ(f->nInstructions(), 5);
 }
 
+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]);
+        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]);
+        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]);
+        CNOT(q[1], q[0]);
+        H(q[0]);
+        H(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");
+    auto circuitCopy = gateRegistry->createComposite("__COPY__");
+    for (size_t i = 0; i < f->nInstructions(); ++i)
+    {
+        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";
+
+    // Validate using gate fusion:
+    auto fuser = xacc::getService<xacc::quantum::GateFuser>("default");
+    fuser->initialize(circuitCopy);
+    const Eigen::MatrixXcd uMatOriginal = fuser->calcFusedGate(5);
+    fuser->initialize(f);
+    const Eigen::MatrixXcd uMatAfter = fuser->calcFusedGate(5);
+    // Compensate any global phase differences.
+    // Find index of the largest element:
+    size_t colIdx = 0;
+    size_t rowIdx = 0;
+    double maxVal = std::abs(uMatAfter(0,0));
+    for (size_t i = 0; i < uMatAfter.rows(); ++i)
+    {
+        for (size_t j = 0; j < uMatAfter.cols(); ++j)
+        {
+            if (std::abs(uMatAfter(i,j)) > maxVal)
+            {
+                maxVal = std::abs(uMatAfter(i,j));
+                colIdx = j;
+                rowIdx = i;
+            }
+        }
+    }
+
+    const std::complex<double> globalFactor = uMatOriginal(rowIdx, colIdx) / uMatAfter(rowIdx, colIdx);
+    auto uMatAfterFixedPhase = globalFactor * uMatAfter;
+    for (size_t i = 0; i < uMatAfter.rows(); ++i)
+    {
+        for (size_t j = 0; j < uMatAfter.cols(); ++j)
+        {
+            EXPECT_NEAR(uMatAfterFixedPhase(i, j).real(), uMatOriginal(i, j).real(), 1e-9);
+            EXPECT_NEAR(uMatAfterFixedPhase(i, j).imag(), uMatOriginal(i, j).imag(), 1e-9);
+        }
+    }
+}
 
 int main(int argc, char **argv) {
   xacc::Initialize(argc, argv);