Mostly completed QPE

- Refactored C-U to become a Circuit class.

- Refactored QPE implementation.

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

quantum/plugins/algorithms/qpe/ControlledGateApplicator.cpp

@@ -13,7 +13,37 @@
 #include "ControlledGateApplicator.hpp"
 
 namespace xacc {
-std::shared_ptr<xacc::CompositeInstruction> ControlledGateApplicator::applyControl(const std::shared_ptr<xacc::CompositeInstruction>& in_program, int in_ctrlIdx)
+namespace circuits {
+
+bool ControlledU::expand(const xacc::HeterogeneousMap& runtimeOptions) 
+{
+    if (!runtimeOptions.keyExists<int>("control-idx")) 
+    {
+        xacc::error("'control-idx' is required.");
+        return false;
+    }
+
+    const auto ctrlIdx = runtimeOptions.get<int>("control-idx");
+    
+    if (!runtimeOptions.pointerLikeExists<CompositeInstruction>("U")) 
+    {
+        xacc::error("'U' composite is required.");
+        return false;
+    }
+    
+    auto uComposite = std::shared_ptr<CompositeInstruction>(
+        runtimeOptions.getPointerLike<CompositeInstruction>("U"), 
+        xacc::empty_delete<CompositeInstruction>());
+    auto ctrlU = applyControl(uComposite, ctrlIdx);
+
+    for (int instId = 0; instId < ctrlU->nInstructions(); ++instId)
+    {
+        addInstruction(ctrlU->getInstruction(instId)->clone());
+    }
+    return true;
+}
+
+std::shared_ptr<xacc::CompositeInstruction> ControlledU::applyControl(const std::shared_ptr<xacc::CompositeInstruction>& in_program, int in_ctrlIdx)
 {
     m_gateProvider = xacc::getService<xacc::IRProvider>("quantum");
     m_composite = m_gateProvider->createComposite("CTRL_" + in_program->name() + "_" + std::to_string(in_ctrlIdx));
@@ -31,7 +61,7 @@ std::shared_ptr<xacc::CompositeInstruction> ControlledGateApplicator::applyContr
     return m_composite;
 }
 
-void ControlledGateApplicator::visit(Hadamard& h) 
+void ControlledU::visit(Hadamard& h) 
 {
     if (h.bits()[0] == m_ctrlIdx)
     {
@@ -45,7 +75,7 @@ void ControlledGateApplicator::visit(Hadamard& h)
     }
 }
 
-void ControlledGateApplicator::visit(X& x) 
+void ControlledU::visit(X& x) 
 {
     if (x.bits()[0] == m_ctrlIdx)
     {
@@ -59,7 +89,7 @@ void ControlledGateApplicator::visit(X& x)
     }
 }
 
-void ControlledGateApplicator::visit(Y& y) 
+void ControlledU::visit(Y& y) 
 {
     if (y.bits()[0] == m_ctrlIdx)
     {
@@ -73,7 +103,7 @@ void ControlledGateApplicator::visit(Y& y)
     }
 }
 
-void ControlledGateApplicator::visit(Z& z) 
+void ControlledU::visit(Z& z) 
 {
     if (z.bits()[0] == m_ctrlIdx)
     {
@@ -87,12 +117,34 @@ void ControlledGateApplicator::visit(Z& z)
     }
 }
 
-void ControlledGateApplicator::visit(CNOT& cnot) 
+void ControlledU::visit(CNOT& cnot) 
 {
-   // TODO: CCNOT
+    // CCNOT gate:
+    // We now has two control bits
+    const auto ctrlIdx1 = cnot.bits()[0];
+    const auto ctrlIdx2 = m_ctrlIdx;
+    // Target qubit
+    const auto targetIdx = cnot.bits()[1];
+    
+    m_composite->addInstruction(m_gateProvider->createInstruction("H", { targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx1, targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("Tdg", { targetIdx }));    
+    m_composite->addInstruction(m_gateProvider->createInstruction("Tdg", { targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx2, targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("T", { targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx1, targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("Tdg", { targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx2, targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("T", { ctrlIdx1 }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("T", { targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("H", { targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx2, ctrlIdx1 }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("T", { ctrlIdx2 }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("Tdg", { ctrlIdx1 }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx2, ctrlIdx1 }));
 }
 
-void ControlledGateApplicator::visit(Rx& rx) 
+void ControlledU::visit(Rx& rx) 
 {
     // CRn(theta) = Rn(theta/2) - CX - Rn(-theta/2) - CX
     if (rx.bits()[0] == m_ctrlIdx)
@@ -110,7 +162,7 @@ void ControlledGateApplicator::visit(Rx& rx)
     }
 }
 
-void ControlledGateApplicator::visit(Ry& ry) 
+void ControlledU::visit(Ry& ry) 
 {
     // CRn(theta) = Rn(theta/2) - CX - Rn(-theta/2) - CX
     if (ry.bits()[0] == m_ctrlIdx)
@@ -129,7 +181,7 @@ void ControlledGateApplicator::visit(Ry& ry)
 }
 
 
-void ControlledGateApplicator::visit(Rz& rz) 
+void ControlledU::visit(Rz& rz) 
 {
     // CRn(theta) = Rn(theta/2) - CX - Rn(-theta/2) - CX
     if (rz.bits()[0] == m_ctrlIdx)
@@ -145,66 +197,69 @@ void ControlledGateApplicator::visit(Rz& rz)
     }
 }
 
-void ControlledGateApplicator::visit(S& s) 
+void ControlledU::visit(S& s) 
 {
-    // S = Rz(pi/2)
-    auto equivGate = m_gateProvider->createInstruction("Rz", { s.bits()[0] }, { M_PI_2 });
-    equivGate->accept(this);
+    // Ctrl-S = CPhase(pi/2)
+    const auto targetIdx = s.bits()[0];
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { m_ctrlIdx, targetIdx }, { M_PI_2 }));
 }
 
-void ControlledGateApplicator::visit(Sdg& sdg) 
+void ControlledU::visit(Sdg& sdg) 
 {
-    // Sdg = Rz(-pi/2)
-    auto equivGate = m_gateProvider->createInstruction("Rz", { sdg.bits()[0] }, { -M_PI_2 });
-    equivGate->accept(this);
+    // Ctrl-Sdg = CPhase(-pi/2)
+    const auto targetIdx = sdg.bits()[0];
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { m_ctrlIdx, targetIdx }, { -M_PI_2 }));
 }
 
-void ControlledGateApplicator::visit(T& t) 
+void ControlledU::visit(T& t) 
 {
-    // T = Rz(pi/4)
-    auto equivGate = m_gateProvider->createInstruction("Rz", { t.bits()[0] }, { M_PI_4 });
-    equivGate->accept(this);
+    // Ctrl-T = CPhase(pi/4)
+    const auto targetIdx = t.bits()[0];
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { m_ctrlIdx, targetIdx }, { M_PI_4 }));
 }
 
-void ControlledGateApplicator::visit(Tdg& tdg) 
+void ControlledU::visit(Tdg& tdg) 
 {
-    // Tdg = Rz(-pi/4)
-    auto equivGate = m_gateProvider->createInstruction("Rz", { tdg.bits()[0] }, { -M_PI_4 });
-    equivGate->accept(this);
+    // Ctrl-Tdg = CPhase(-pi/4)
+    const auto targetIdx = tdg.bits()[0];
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { m_ctrlIdx, targetIdx }, { -M_PI_4 }));
 }
 
-void ControlledGateApplicator::visit(U& u) 
+void ControlledU::visit(Swap& s) 
 {
-    // TODO
+    CNOT c1(s.bits()), c2(s.bits()[1],s.bits()[0]), c3(s.bits());
+    visit(c1);
+    visit(c2);
+    visit(c3);
 }
 
-void ControlledGateApplicator::visit(CY& cy) 
+void ControlledU::visit(U& u) 
 {
-    // TODO
+    xacc::error("Unsupported!");
 }
 
-void ControlledGateApplicator::visit(CZ& cz) 
+void ControlledU::visit(CY& cy) 
 {
-    // TODO
+    xacc::error("Unsupported!");
 }
 
-void ControlledGateApplicator::visit(CRZ& crz) 
+void ControlledU::visit(CZ& cz) 
 {
-    // TODO
+    xacc::error("Unsupported!");
 }
 
-void ControlledGateApplicator::visit(CH& ch) 
+void ControlledU::visit(CRZ& crz) 
 {
-    // TODO
+    xacc::error("Unsupported!");
 }
 
-void ControlledGateApplicator::visit(Swap& s) 
+void ControlledU::visit(CH& ch) 
 {
-   // TODO
+    xacc::error("Unsupported!");
 }
 
-void ControlledGateApplicator::visit(CPhase& cphase) 
+void ControlledU::visit(CPhase& cphase) 
 {
-   // TODO
+    xacc::error("Unsupported!");
 }
-}
\ No newline at end of file
+}}
\ No newline at end of file

quantum/plugins/algorithms/qpe/ControlledGateApplicator.hpp

@@ -16,12 +16,18 @@
 using namespace xacc::quantum;
 
 namespace xacc {
-class ControlledGateApplicator: public AllGateVisitor
+namespace circuits {
+
+class ControlledU: public AllGateVisitor, public quantum::Circuit
 {
 public:
-    // Apply *single* control on the input composite.
-    // Returns the new composite.
-    std::shared_ptr<xacc::CompositeInstruction> applyControl(const std::shared_ptr<xacc::CompositeInstruction>& in_program, int in_ctrlIdx);
+    ControlledU() : Circuit("C-U") {}
+    bool expand(const xacc::HeterogeneousMap& runtimeOptions) override;
+    // Input: The composite "U" and the control Idx.
+    // Control Idx must *not* be one of the qubits that U is acting on.
+    const std::vector<std::string> requiredKeys() override { return { "U", "control-idx" }; }
+    DEFINE_CLONE(ControlledU);
+
     // AllGateVisitor implementation
     void visit(Hadamard& h) override;
     void visit(CNOT& cnot) override;
@@ -49,9 +55,15 @@ public:
     void visit(IfStmt& ifStmt) override { xacc::error("Unsupported!"); }
     void visit(fSim& fsim) override { xacc::error("Unsupported!"); }
     void visit(iSwap& isw) override { xacc::error("Unsupported!"); }
+
+private:
+    // Apply *single* control on the input composite.
+    // Returns the new composite.
+    std::shared_ptr<xacc::CompositeInstruction> applyControl(const std::shared_ptr<xacc::CompositeInstruction>& in_program, int in_ctrlIdx);
+
 private:
     std::shared_ptr<xacc::CompositeInstruction> m_composite;
     std::shared_ptr<xacc::IRProvider> m_gateProvider;
     size_t m_ctrlIdx;
 };
-}
\ No newline at end of file
+}}
\ No newline at end of file

quantum/plugins/algorithms/qpe/QpeActivator.cpp

@@ -11,6 +11,7 @@
  *   Thien Nguyen - initial API and implementation
  *******************************************************************************/
 #include "qpe.hpp"
+#include "ControlledGateApplicator.hpp"
 
 #include "cppmicroservices/BundleActivator.h"
 #include "cppmicroservices/BundleContext.h"
@@ -35,6 +36,7 @@ public:
   void Start(BundleContext context) {
     auto c = std::make_shared<xacc::algorithm::QuantumPhaseEstimation>();
     context.RegisterService<xacc::Algorithm>(c);
+    context.RegisterService<xacc::Instruction>(std::make_shared<xacc::circuits::ControlledU>());
   }
 
   /**

quantum/plugins/algorithms/qpe/qpe.cpp

@@ -17,7 +17,6 @@
 #include "Circuit.hpp"
 #include <cassert>
 #include <iomanip>
-#include "ControlledGateApplicator.hpp" 
 
 namespace xacc {
 namespace algorithm {
@@ -54,78 +53,104 @@ const std::vector<std::string> QuantumPhaseEstimation::requiredParameters() cons
 
 void QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const 
 {
-    if (buffer->size() < 1)
+    // Calculate the number of qubits that are used by the oracle:
+    const size_t nbOracleBits = m_oracle->uniqueBits().size();
+
+    if (nbOracleBits < 1)
+    {
+        std::cout << "Invalid oracle Composite.\n";
+    }
+    
+    // Buffer must contain additional qubits for phase est. result
+    if (buffer->size() <= nbOracleBits)
     {
-        std::cout << "Buffer must have more than 1 qubits.\n";
+        std::cout << "Buffer must have more than " << nbOracleBits << " qubits.\n";
     }
     
     auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
     auto qpeKernel = gateRegistry->createComposite("QuantumPhaseEstKernel");
 
     // Bit precision: the number of extra qubits that were allocated in the input buffer.
-    const auto bitPrecision = buffer->size() - 1;
+    const auto bitPrecision = buffer->size() - nbOracleBits;
     std::cout << "Phase estimation precision = " << bitPrecision << " bits.\n";
-    // Convention: q[0] is the main qubit
-    // q[1]..q[n] are the phase result qubits
+    
+    // Map the oracle/state-prep circuit to the qubit range after the phase estimation register:
+    // q[0] -> q[bitPrecision]
+    // q[1] -> q[bitPrecision + 1], etc...
+    const auto mapPrimaryQubits = [&bitPrecision](const std::shared_ptr<xacc::CompositeInstruction>& in_composite){
+        InstructionIterator it(in_composite);
+        while (it.hasNext())
+        {
+            auto nextInst = it.next();
+            if (nextInst->isEnabled())
+            {
+                auto currentBits = nextInst->bits();
+                for (auto& bit: currentBits)
+                {
+                    bit = bit + bitPrecision;
+                }
+                nextInst->setBits(currentBits);
+            }
+        }
+    };
+    
+    // Convention: q[0] - q[precision - 1] is phase-estimation result register.
     // Hadamard on all ancilla/result qubits
-    for (size_t i = 1; i < buffer->size(); ++i)
+    for (size_t i = 0; i < bitPrecision; ++i)
     {
         qpeKernel->addInstruction(gateRegistry->createInstruction("H", { i }));
     }
 
-    // Prepare q[0] in the eigenstate:
+    // Prepare the eigenstate:
     // Note: user must provide 'state-preparation' composite 
     // to transform |0> state to the eigenstate.
     // Otherwise, assume |0> is the eigenstate.
-    if (!m_params.pointerLikeExists<CompositeInstruction>("state-preparation")) 
+    if (m_params.pointerLikeExists<CompositeInstruction>("state-preparation")) 
     {
-        auto statePrep = m_params.getPointerLike<CompositeInstruction>("state-preparation");
-        if (statePrep->uniqueBits().size() != 1)
+        auto statePrep = m_params.getPointerLike<CompositeInstruction>("state-preparation");        
+        mapPrimaryQubits(std::shared_ptr<CompositeInstruction>(statePrep, xacc::empty_delete<CompositeInstruction>()));
+
+        if (statePrep->uniqueBits().size() > nbOracleBits)
         {
-            xacc::error("'state-preparation' circuit should only contain one qubit.");
+            xacc::error("'state-preparation' circuit cannot operate on more qubits than the oracle.");
             return;
         }
+
         // Add state preparation composite
         qpeKernel->addInstructions(statePrep->getInstructions());
     }
 
+    auto oracleSharedPtr = std::shared_ptr<CompositeInstruction>(m_oracle, xacc::empty_delete<CompositeInstruction>());
+    mapPrimaryQubits(oracleSharedPtr);
+    
     // Controlled-oracle application
-    ControlledGateApplicator gateApplicator;
-    for (size_t i = 1; i < buffer->size(); ++i)
+    for (size_t i = 0; i < bitPrecision; ++i)
     {
-        // q1: U; q2: U^2; q3: U^4; etc.
-        const int nbCalls = 1 << (i - 1);
-        
-        auto ctrlKernel = gateApplicator.applyControl(std::shared_ptr<CompositeInstruction>(m_oracle, xacc::empty_delete<CompositeInstruction>()), i);
+        // q0: U; q1: U^2; q2: U^4; etc.
+        const int nbCalls = 1 << i;
+        auto ctrlKernel = std::dynamic_pointer_cast<CompositeInstruction>(xacc::getService<Instruction>("C-U"));
+        ctrlKernel->expand( { 
+            std::make_pair("U",  m_oracle),
+            std::make_pair("control-idx",  static_cast<int>(i)),
+        });            
         // Apply C-U^n
         for (int count = 0; count < nbCalls; count++)
         {
-            qpeKernel->addInstructions(ctrlKernel->getInstructions());
-        }
-    }
-    
-    // IQFT on q[1]-q[n]
-    auto iqft = std::dynamic_pointer_cast<CompositeInstruction>(xacc::getService<Instruction>("iqft"));
-    iqft->expand( { std::make_pair("nq", bitPrecision) });
-    // We need to shift the qubit index up by 1 
-    InstructionIterator it(iqft);
-    while (it.hasNext())
-    {
-        auto nextInst = it.next();
-        if (nextInst->isEnabled())
-        {
-            auto currentBits = nextInst->bits();
-            for (auto& bit: currentBits)
+            for (int instId = 0; instId < ctrlKernel->nInstructions(); ++instId)
             {
-                bit = bit + 1;
+                // We need to clone the instruction since it'll be repeated.
+                qpeKernel->addInstruction(ctrlKernel->getInstruction(instId)->clone());
             }
-            nextInst->setBits(currentBits);
         }
     }
+    
+    // IQFT on the phase estimation register.
+    auto iqft = std::dynamic_pointer_cast<CompositeInstruction>(xacc::getService<Instruction>("iqft"));
+    iqft->expand( { std::make_pair("nq", static_cast<int>(bitPrecision)) });    
     qpeKernel->addInstructions(iqft->getInstructions());
-
+    
     // Measure the ancilla/result qubits
-    for (size_t i = 1; i < buffer->size(); ++i)
+    for (size_t i = 0; i < bitPrecision; ++i)
     {
         qpeKernel->addInstruction(gateRegistry->createInstruction("Measure", { i }));
     }
@@ -133,7 +158,7 @@ void QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> bu
     // DEBUG: 
     std::cout << "QPE kernel:\n" << qpeKernel->toString() << "\n\n";
 
-    // TODO: execute
+    m_qpu->execute(buffer, qpeKernel);
 }
 
 std::vector<double> QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<double>& x) 

quantum/plugins/algorithms/qpe/tests/QpeTester.cpp

@@ -20,14 +20,16 @@ using namespace xacc;
 
 TEST(QpeTester, checkSimple) 
 {
-  auto acc = xacc::getAccelerator("qpp");
+  auto acc = xacc::getAccelerator("qpp", {std::make_pair("shots", 1024)});
   // Test case: T gate, eigenstate = |1>
   // 3-bit precision
   auto buffer = xacc::qalloc(4);
   auto qpe = xacc::getService<Algorithm>("QPE");
   auto compiler = xacc::getCompiler("xasm");
   
-  // Oracle = T gate
+  // Oracle = T gate 
+  // |1> => exp(i*pi/4) |1>
+  // Expected result = |100> = |4>
   auto oracle = compiler->compile(R"(__qpu__ void oracle(qbit q) {
     T(q[0]); 
   })", nullptr)->getComposite("oracle");
@@ -43,6 +45,43 @@ TEST(QpeTester, checkSimple)
                     std::make_pair("state-preparation", statePrep)
                   }));
   qpe->execute(buffer);
+
+  const auto result = buffer->computeMeasurementProbability("100");
+  EXPECT_NEAR(result, 1.0, 0.01);
+}
+
+TEST(QpeTester, checkMultiQubitOracle) 
+{
+  auto acc = xacc::getAccelerator("qpp", {std::make_pair("shots", 1024)});
+  // Test case: T gate on both qubits, eigenstate = |11>
+  // i.e. TT |11> = exp(i*pi/4)* exp(i*pi/4)*|11> = exp(i*pi/2)*|11>
+  // => the expected answer is 2 = 010
+
+  // 3-bit precision => 5 qubits in total
+  auto buffer = xacc::qalloc(5);
+  auto qpe = xacc::getService<Algorithm>("QPE");
+  auto compiler = xacc::getCompiler("xasm");
+  
+  // Oracle: T gate on both qubits
+  auto oracle = compiler->compile(R"(__qpu__ void oracle(qbit q) {
+    T(q[0]); 
+    T(q[1]); 
+  })", nullptr)->getComposite("oracle");
+
+  // Eigenstate preparation = |11> state
+  auto statePrep = compiler->compile(R"(__qpu__ void prep1(qbit q) {
+    X(q[0]); 
+    X(q[1]); 
+  })", nullptr)->getComposite("prep1");  
+  
+  EXPECT_TRUE(qpe->initialize({
+                    std::make_pair("accelerator", acc),
+                    std::make_pair("oracle", oracle),
+                    std::make_pair("state-preparation", statePrep)
+                  }));
+  qpe->execute(buffer);
+  const auto result = buffer->computeMeasurementProbability("010");
+  EXPECT_NEAR(result, 1.0, 0.01);
 }
 
 int main(int argc, char **argv)