Completed QPE

- Implemeted the remaining CU gates.

- Added Cpp and Python examples.
Signed-off-by: Thien Nguyen <nguyentm@ornl.gov>

python/examples/quantum_phase_estimation.py

+import xacc,sys, numpy as np
+
+# Get access to the desired QPU and
+# allocate some qubits to run on
+qpu = xacc.getAccelerator('qpp',  { 'shots': 4096 })
+
+# In this example: we want to estimate the *phase* of an arbitrary 'oracle'
+# i.e. Oracle(|State>) = exp(i*Phase)*|State>
+# and we need to estimate that Phase.
+
+# The oracle is a simple T gate, and the eigenstate is |1>
+# T|1> = e^(i*pi/4)|1> 
+# The phase value of pi/4 = 2pi * (1/8)
+# i.e. if we use a 3-bit register for estimation, 
+# we will get the correct answer of 1 deterministically.
+
+xacc.qasm('''.compiler xasm
+.circuit oracle
+.qbit q
+T(q[0]);
+''')
+oracle = xacc.getCompiled('oracle')
+
+# We need to prepare the eigenstate |1>
+xacc.qasm('''.compiler xasm
+.circuit prep
+.qbit q
+X(q[0]);
+''')
+statePrep = xacc.getCompiled('prep')
+
+# We need 4 qubits (3-bit precision)
+buffer = xacc.qalloc(4)
+
+# Create the QPE algorithm
+qpe = xacc.getAlgorithm('QPE', {
+                        'accelerator': qpu,
+                        'oracle': oracle,
+                        'state-preparation': statePrep
+                        })
+
+qpe.execute(buffer)
+# We should only get the bit string of |100> = 1
+# i.e. phase value of 1/2^3 = 1/8.
+print(buffer)
\ No newline at end of file

quantum/examples/qasm/CMakeLists.txt

@@ -35,3 +35,6 @@ target_link_libraries(nah_ucc3 PRIVATE xacc xacc-quantum-gate)
 
 add_executable(qaoa_qubo solve_qubo_with_qaoa.cpp)
 target_link_libraries(qaoa_qubo PRIVATE xacc xacc-quantum-gate)
+
+add_executable(quantum_phase_estimation quantum_phase_estimation.cpp)
+target_link_libraries(quantum_phase_estimation PRIVATE xacc xacc-quantum-gate)

quantum/examples/qasm/quantum_phase_estimation.cpp

+/*******************************************************************************
+ * Copyright (c) 2020 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 "xacc.hpp"
+#include "xacc_observable.hpp"
+#include "xacc_service.hpp"
+#include <random>
+
+int main(int argc, char **argv) {
+  xacc::Initialize(argc, argv);
+  // Accelerator:
+  auto acc = xacc::getAccelerator("qpp", {std::make_pair("shots", 4096)});
+  
+  // In this example: we want to estimate the *phase* of an arbitrary 'oracle'
+  // i.e. Oracle(|State>) = exp(i*Phase)*|State>
+  // and we need to estimate that Phase.
+
+  // Oracle: CPhase(theta) or CU1(theta) which is defined as
+  // 1 0 0 0
+  // 0 1 0 0
+  // 0 0 1 0
+  // 0 0 0 e^(i*theta)
+  // The eigenstate is |11>; i.e. CPhase(theta)|11> = e^(i*theta)|11>
+
+  // Since this oracle operates on 2 qubits, we need to add more qubits to the buffer.
+  // The more qubits we have, the more accurate the estimate.
+  // Resolution := 2^(number qubits in the calculation register).
+  // 5-bit precision => 7 qubits in total
+  auto buffer = xacc::qalloc(7);
+  auto qpe = xacc::getService<xacc::Algorithm>("QPE");
+  auto compiler = xacc::getCompiler("xasm");
+  
+  // Create oracle: CPhase gate with theta = 2pi/3
+  // i.e. the phase value to estimate is 1/3 ~ 0.33333.
+  auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
+  auto oracle = gateRegistry->createComposite("oracle");
+  oracle->addInstruction(gateRegistry->createInstruction("CPhase", { 0, 1 }, { 2.0 * M_PI/ 3.0 }));
+
+  // Eigenstate preparation = |11> state
+  auto statePrep = compiler->compile(R"(__qpu__ void prep1(qbit q) {
+    X(q[0]); 
+    X(q[1]); 
+  })", nullptr)->getComposite("prep1");  
+  
+  // Initialize the Quantum Phase Estimation:
+  qpe->initialize({
+                    std::make_pair("accelerator", acc),
+                    std::make_pair("oracle", oracle),
+                    std::make_pair("state-preparation", statePrep)
+                  });
+  
+  // Run the algorithm
+  qpe->execute(buffer);
+  // Expected result: 
+  // The factor here is 2^5 (precision) = 32
+  // we expect the two most-likely bitstring is 10 and 11
+  // i.e. the true result is between 10/32 = 0.3125 and 11/32 = 0.34375  
+  std::cout << "Probability of the two most-likely bitstrings 10 (theta = 0.3125) and 11 (theta = 0.34375 ): \n";
+  std::cout << "Probability of |11010> (11) = " << buffer->computeMeasurementProbability("11010") << "\n";
+  std::cout << "Probability of |01010> (10) = " << buffer->computeMeasurementProbability("01010") << "\n";
+
+  xacc::Finalize();
+}
\ No newline at end of file

quantum/plugins/algorithms/qpe/ControlledGateApplicator.cpp

@@ -120,7 +120,7 @@ void ControlledU::visit(Z& z)
 void ControlledU::visit(CNOT& cnot) 
 {
     // CCNOT gate:
-    // We now has two control bits
+    // We now have two control bits
     const auto ctrlIdx1 = cnot.bits()[0];
     const auto ctrlIdx2 = m_ctrlIdx;
     // Target qubit
@@ -183,7 +183,6 @@ void ControlledU::visit(Ry& ry)
 
 void ControlledU::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).");
@@ -227,6 +226,9 @@ void ControlledU::visit(Tdg& tdg)
 
 void ControlledU::visit(Swap& s) 
 {
+    // Fredkin gate: controlled-swap
+    // Decompose Swap gate into CNOT gates;
+    // then re-visit those CNOT gates (becoming CCNOT).
     CNOT c1(s.bits()), c2(s.bits()[1],s.bits()[0]), c3(s.bits());
     visit(c1);
     visit(c2);
@@ -235,31 +237,93 @@ void ControlledU::visit(Swap& s)
 
 void ControlledU::visit(U& u) 
 {
-    xacc::error("Unsupported!");
+    // U(theta, phi, lambda) := Rz(phi)Ry(theta)Rz(lambda)
+    const auto theta = InstructionParameterToDouble(u.getParameter(0));
+    const auto phi = InstructionParameterToDouble(u.getParameter(1));
+    const auto lambda = InstructionParameterToDouble(u.getParameter(2));
+    Ry ry1(u.bits()[0], theta);
+    Rz rz1(u.bits()[0], lambda);
+    Rz rz2(u.bits()[0], phi);
+    // Revisit these gates: Rn -> CRn
+    visit(rz1);
+    visit(ry1);
+    visit(rz2);
 }
 
 void ControlledU::visit(CY& cy) 
 {
-    xacc::error("Unsupported!");
+    // controlled-Y = Sdg(target) - CX - S(target)
+    CNOT c1(cy.bits());
+    Sdg sdg(cy.bits()[1]);
+    S s(cy.bits()[1]);
+    // Revisit these gates: CNOT->CCNOT; S -> CPhase
+    visit(sdg);
+    visit(c1);
+    visit(s);
 }
 
 void ControlledU::visit(CZ& cz) 
 {
-    xacc::error("Unsupported!");
+    // CZ = H(target) - CX - H(target) 
+    CNOT c1(cz.bits());
+    Hadamard h1(cz.bits()[1]);
+    Hadamard h2(cz.bits()[1]);
+
+    // Revisit these gates: CNOT->CCNOT; H -> CH
+    visit(h1);
+    visit(c1);
+    visit(h2);
 }
 
 void ControlledU::visit(CRZ& crz) 
 {
-    xacc::error("Unsupported!");
+    const auto theta = InstructionParameterToDouble(crz.getParameter(0));
+    // Decompose
+    Rz rz1(crz.bits()[1], theta/2);
+    CNOT c1(crz.bits());
+    Rz rz2(crz.bits()[1], -theta/2);
+    CNOT c2(crz.bits());
+    // Revisit:
+    visit(rz1);
+    visit(c1);
+    visit(rz2);
+    visit(c2);
 }
 
 void ControlledU::visit(CH& ch) 
 {
-    xacc::error("Unsupported!");
+    // controlled-H = Ry(pi/4, target) - CX - Ry(-pi/4, target)
+    CNOT c1(ch.bits());
+    Ry ry1(ch.bits()[1], M_PI_4);
+    Ry ry2(ch.bits()[1], -M_PI_4);
+    visit(ry1);
+    visit(c1);
+    visit(ry2);
 }
 
 void ControlledU::visit(CPhase& cphase) 
 {
-    xacc::error("Unsupported!");
+    // This is effectively CRz but due to the global phase difference
+    // in the matrix definitions of U1 and Rz, 
+    // CU1 (i.e. CPhase) and CRz are different gates with a relative phase difference.
+    // Ref:
+    // ccu1(t, ctrl1, ctrl2, target) =
+    // cu1(t/2, ctrl1, ctrl2)
+    // cx(ctrl2, target)
+    // cu1(-t/2, ctrl1, target)
+    // cx(ctrl2, target)
+    // cu1(t/2, ctrl1, target)
+
+    const auto ctrlIdx1 = cphase.bits()[0];
+    const auto ctrlIdx2 = m_ctrlIdx;
+    // Target qubit
+    const auto targetIdx = cphase.bits()[1];
+    // Angle
+    const auto angle = InstructionParameterToDouble(cphase.getParameter(0));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { ctrlIdx1, ctrlIdx2 }, { angle/2 }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx2, targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { ctrlIdx1, targetIdx }, { -angle/2 }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CX", { ctrlIdx2, targetIdx }));
+    m_composite->addInstruction(m_gateProvider->createInstruction("CPhase", { ctrlIdx1, targetIdx }, { angle/2 }));
 }
 }}
\ No newline at end of file

quantum/plugins/algorithms/qpe/qpe.cpp

@@ -155,15 +155,14 @@ void QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> bu
         qpeKernel->addInstruction(gateRegistry->createInstruction("Measure", { i }));
     }
     
-    // DEBUG: 
-    std::cout << "QPE kernel:\n" << qpeKernel->toString() << "\n\n";
-
+    buffer->addExtraInfo("qpe-kernel", qpeKernel->toString());
     m_qpu->execute(buffer, qpeKernel);
 }
 
 std::vector<double> QuantumPhaseEstimation::execute(const std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<double>& x) 
 {
-    // TODO
+    // We don't have this for QPE.
+    xacc::error("This method is unsupported!");
     return { };
 }
 } // namespace algorithm