Work on QPE algorithm

- Template for applying controlled oracle.

- Contruct the basic circuit for the algorithm.
Signed-off-by: Thien Nguyen <nguyentm@ornl.gov>

quantum/plugins/algorithms/qpe/ControlledGateApplicator.cpp

+/*******************************************************************************
+ * Copyright (c) 2019 UT-Battelle, LLC.
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v1.0 which accompanies this
+ * distribution. The Eclipse Public License is available at
+ * http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution
+ *License is available at https://eclipse.org/org/documents/edl-v10.php
+ *
+ * Contributors:
+ *   Thien Nguyen - initial API and implementation
+ *******************************************************************************/
+#include "ControlledGateApplicator.hpp"
+
+namespace xacc {
+std::shared_ptr<xacc::CompositeInstruction> ControlledGateApplicator::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));
+    m_ctrlIdx = in_ctrlIdx;
+    InstructionIterator it(in_program);
+    while (it.hasNext())
+    {
+        auto nextInst = it.next();
+        if (nextInst->isEnabled())
+        {
+            nextInst->accept(this);
+        }
+    }
+
+    return m_composite;
+}
+
+void ControlledGateApplicator::visit(Hadamard& h) 
+{
+    if (h.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = h.bits()[0];
+        // CH
+        m_composite->addInstruction(m_gateProvider->createInstruction("CH", { m_ctrlIdx, targetIdx }));
+    }
+}
+
+void ControlledGateApplicator::visit(X& x) 
+{
+    if (x.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = x.bits()[0];
+        // CNOT
+        m_composite->addInstruction(m_gateProvider->createInstruction("CX", { m_ctrlIdx, targetIdx }));
+    }
+}
+
+void ControlledGateApplicator::visit(Y& y) 
+{
+    if (y.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = y.bits()[0];
+        // CY
+        m_composite->addInstruction(m_gateProvider->createInstruction("CY", { m_ctrlIdx, targetIdx }));
+    }
+}
+
+void ControlledGateApplicator::visit(Z& z) 
+{
+    if (z.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = z.bits()[0];
+        // CZ
+        m_composite->addInstruction(m_gateProvider->createInstruction("CZ", { m_ctrlIdx, targetIdx }));
+    }
+}
+
+void ControlledGateApplicator::visit(CNOT& cnot) 
+{
+   // TODO: CCNOT
+}
+
+void ControlledGateApplicator::visit(Rx& rx) 
+{
+    // CRn(theta) = Rn(theta/2) - CX - Rn(-theta/2) - CX
+    if (rx.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = rx.bits()[0];
+        const auto angle = InstructionParameterToDouble(rx.getParameter(0));
+        m_composite->addInstruction(m_gateProvider->createInstruction("Rx", { targetIdx }, { angle/ 2.0 }));
+        m_composite->addInstruction(m_gateProvider->createInstruction("CX", { m_ctrlIdx, targetIdx }));
+        m_composite->addInstruction(m_gateProvider->createInstruction("Rx", { targetIdx }, { -angle/ 2.0 }));
+        m_composite->addInstruction(m_gateProvider->createInstruction("CX", { m_ctrlIdx, targetIdx }));
+    }
+}
+
+void ControlledGateApplicator::visit(Ry& ry) 
+{
+    // CRn(theta) = Rn(theta/2) - CX - Rn(-theta/2) - CX
+    if (ry.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = ry.bits()[0];
+        const auto angle = InstructionParameterToDouble(ry.getParameter(0));
+        m_composite->addInstruction(m_gateProvider->createInstruction("Ry", { targetIdx }, { angle/ 2.0 }));
+        m_composite->addInstruction(m_gateProvider->createInstruction("CX", { m_ctrlIdx, targetIdx }));
+        m_composite->addInstruction(m_gateProvider->createInstruction("Ry", { targetIdx }, { -angle/ 2.0 }));
+        m_composite->addInstruction(m_gateProvider->createInstruction("CX", { m_ctrlIdx, targetIdx }));
+    }
+}
+
+
+void ControlledGateApplicator::visit(Rz& rz) 
+{
+    // CRn(theta) = Rn(theta/2) - CX - Rn(-theta/2) - CX
+    if (rz.bits()[0] == m_ctrlIdx)
+    {
+        xacc::error("Control bit must be different from target qubit(s).");
+    }
+    else
+    {
+        const auto targetIdx = rz.bits()[0];
+        const auto angle = InstructionParameterToDouble(rz.getParameter(0));
+        // CRz
+        m_composite->addInstruction(m_gateProvider->createInstruction("CRZ", { m_ctrlIdx, targetIdx }, { angle }));
+    }
+}
+
+void ControlledGateApplicator::visit(S& s) 
+{
+    // S = Rz(pi/2)
+    auto equivGate = m_gateProvider->createInstruction("Rz", { s.bits()[0] }, { M_PI_2 });
+    equivGate->accept(this);
+}
+
+void ControlledGateApplicator::visit(Sdg& sdg) 
+{
+    // Sdg = Rz(-pi/2)
+    auto equivGate = m_gateProvider->createInstruction("Rz", { sdg.bits()[0] }, { -M_PI_2 });
+    equivGate->accept(this);
+}
+
+void ControlledGateApplicator::visit(T& t) 
+{
+    // T = Rz(pi/4)
+    auto equivGate = m_gateProvider->createInstruction("Rz", { t.bits()[0] }, { M_PI_4 });
+    equivGate->accept(this);
+}
+
+void ControlledGateApplicator::visit(Tdg& tdg) 
+{
+    // Tdg = Rz(-pi/4)
+    auto equivGate = m_gateProvider->createInstruction("Rz", { tdg.bits()[0] }, { -M_PI_4 });
+    equivGate->accept(this);
+}
+
+void ControlledGateApplicator::visit(U& u) 
+{
+    // TODO
+}
+
+void ControlledGateApplicator::visit(CY& cy) 
+{
+    // TODO
+}
+
+void ControlledGateApplicator::visit(CZ& cz) 
+{
+    // TODO
+}
+
+void ControlledGateApplicator::visit(CRZ& crz) 
+{
+    // TODO
+}
+
+void ControlledGateApplicator::visit(CH& ch) 
+{
+    // TODO
+}
+
+void ControlledGateApplicator::visit(Swap& s) 
+{
+   // TODO
+}
+
+void ControlledGateApplicator::visit(CPhase& cphase) 
+{
+   // TODO
+}
+}
\ No newline at end of file

quantum/plugins/algorithms/qpe/ControlledGateApplicator.hpp

+/*******************************************************************************
+ * Copyright (c) 2019 UT-Battelle, LLC.
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v1.0 which accompanies this
+ * distribution. The Eclipse Public License is available at
+ * http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution
+ *License is available at https://eclipse.org/org/documents/edl-v10.php
+ *
+ * Contributors:
+ *   Thien Nguyen - initial API and implementation
+ *******************************************************************************/
+#pragma once
+
+#include "AllGateVisitor.hpp"
+using namespace xacc::quantum;
+
+namespace xacc {
+class ControlledGateApplicator: public AllGateVisitor
+{
+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);
+    // AllGateVisitor implementation
+    void visit(Hadamard& h) override;
+    void visit(CNOT& cnot) override;
+    void visit(Rz& rz) override;
+    void visit(Ry& ry) override;
+    void visit(Rx& rx) override;
+    void visit(X& x) override;
+    void visit(Y& y) override;
+    void visit(Z& z) override;
+    void visit(CY& cy) override;
+    void visit(CZ& cz) override;
+    void visit(Swap& s) override;
+    void visit(CRZ& crz) override;
+    void visit(CH& ch) override;
+    void visit(S& s) override;
+    void visit(Sdg& sdg) override;
+    void visit(T& t) override;
+    void visit(Tdg& tdg) override;
+    void visit(CPhase& cphase) override;
+    void visit(U& u) override;
+    // Nothing to do
+    void visit(Identity& i) override {}
+    // These are unsupported.
+    void visit(Measure& measure) override { xacc::error("Unsupported!"); }
+    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:
+    std::shared_ptr<xacc::CompositeInstruction> m_composite;
+    std::shared_ptr<xacc::IRProvider> m_gateProvider;
+    size_t m_ctrlIdx;
+};
+}
\ No newline at end of file

quantum/plugins/algorithms/qpe/qpe.cpp

@@ -17,6 +17,7 @@
 #include "Circuit.hpp"
 #include <cassert>
 #include <iomanip>
+#include "ControlledGateApplicator.hpp" 
 
 namespace xacc {
 namespace algorithm {
@@ -58,11 +59,81 @@ void QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> bu
         std::cout << "Buffer must have more than 1 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;
+    std::cout << "Phase estimation precision = " << bitPrecision << " bits.\n";
+    // Convention: q[0] is the main qubit
+    // q[1]..q[n] are the phase result qubits
+    // Hadamard on all ancilla/result qubits
+    for (size_t i = 1; i < buffer->size(); ++i)
+    {
+        qpeKernel->addInstruction(gateRegistry->createInstruction("H", { i }));
+    }
 
+    // Prepare q[0] in 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")) 
+    {
+        auto statePrep = m_params.getPointerLike<CompositeInstruction>("state-preparation");
+        if (statePrep->uniqueBits().size() != 1)
+        {
+            xacc::error("'state-preparation' circuit should only contain one qubit.");
+            return;
+        }
+        // Add state preparation composite
+        qpeKernel->addInstructions(statePrep->getInstructions());
+    }
 
+    // Controlled-oracle application
+    ControlledGateApplicator gateApplicator;
+    for (size_t i = 1; i < buffer->size(); ++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);
+        // 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)
+            {
+                bit = bit + 1;
+            }
+            nextInst->setBits(currentBits);
+        }
+    }
+    qpeKernel->addInstructions(iqft->getInstructions());
 
+    // Measure the ancilla/result qubits
+    for (size_t i = 1; i < buffer->size(); ++i)
+    {
+        qpeKernel->addInstruction(gateRegistry->createInstruction("Measure", { i }));
+    }
+    
+    // DEBUG: 
+    std::cout << "QPE kernel:\n" << qpeKernel->toString() << "\n\n";
+
+    // TODO: execute
 }
 
 std::vector<double> QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<double>& x) 
@@ -70,6 +141,5 @@ std::vector<double> QuantumPhaseEstimation::execute(const std::shared_ptr<Accele
     // TODO
     return { };
 }
-
 } // namespace algorithm
 } // namespace xacc
\ No newline at end of file

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

@@ -14,15 +14,35 @@
 
 #include "xacc.hpp"
 #include "xacc_service.hpp"
-#include "Observable.hpp"
 #include "Algorithm.hpp"
-#include "xacc_observable.hpp"
 
 using namespace xacc;
 
 TEST(QpeTester, checkSimple) 
 {
   auto acc = xacc::getAccelerator("qpp");
+  // 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
+  auto oracle = compiler->compile(R"(__qpu__ void oracle(qbit q) {
+    T(q[0]); 
+  })", nullptr)->getComposite("oracle");
+  
+  // Eigenstate preparation = X gate (|1> state)
+  auto statePrep = compiler->compile(R"(__qpu__ void prep(qbit q) {
+    X(q[0]); 
+  })", nullptr)->getComposite("prep");  
+  
+  EXPECT_TRUE(qpe->initialize({
+                    std::make_pair("accelerator", acc),
+                    std::make_pair("oracle", oracle),
+                    std::make_pair("state-preparation", statePrep)
+                  }));
+  qpe->execute(buffer);
 }
 
 int main(int argc, char **argv)