Implemented bitstring direct sampling for the Qpp Acc

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

quantum/plugins/qpp/accelerator/QppAccelerator.cpp

@@ -11,6 +11,169 @@
  *   Thien Nguyen - initial API and implementation
  *******************************************************************************/
 #include "QppAccelerator.hpp"
+#include <thread>
+#include <mutex>
+
+namespace {
+    inline bool isMeasureGate(const xacc::InstPtr& in_instr)
+    {
+        return (in_instr->name() == "Measure");
+    }
+
+    inline size_t getNumberOfThreads()
+    {
+        static const size_t NB_THREADS = std::thread::hardware_concurrency();
+        return NB_THREADS;
+    }
+
+    inline double generateRandomProbability() 
+    {
+        auto randFunc = std::bind(std::uniform_real_distribution<double>(0, 1), std::mt19937(std::chrono::high_resolution_clock::now().time_since_epoch().count()));
+        return randFunc();
+    }
+
+    bool applyMeasureOp(KetVectorType& io_psi, size_t in_qubitIndex)
+    {
+        const auto N = io_psi.size();
+        const auto k_range = 1ULL << in_qubitIndex;        
+        const double randProbPick = generateRandomProbability();
+        double cumulativeProb = 0.0;
+        size_t stateSelect = 0;
+        //select a state based on cumulative distribution
+        while (cumulativeProb < randProbPick && stateSelect < N)
+        {
+            cumulativeProb += std::norm(io_psi[stateSelect++]);
+        }
+        stateSelect--;
+
+        //take the value of the measured bit
+        bool result = ((stateSelect >> in_qubitIndex) & 1);
+        // Collapse the state vector according to the measurement result
+        double measProb = 0.0;
+        for (size_t g = 0; g < N; g += (k_range * 2))
+        {
+            for (size_t i = 0; i < k_range; ++i)
+            {
+                if ((((i + g) >> in_qubitIndex) & 1) == result)
+                {
+                    measProb += std::norm(io_psi[i + g]);
+                    io_psi[i + g + k_range] = 0.0; 
+                }
+                else
+                {
+                    measProb += std::norm(io_psi[i + g + k_range]);
+                    io_psi[i + g] = 0.0;
+                }
+            }
+        }
+        // Renormalize the state
+        measProb = std::sqrt(measProb);
+        for (size_t g = 0; g < N; g += (k_range * 2))
+        {
+            for (size_t i = 0; i < k_range; i += 1)
+            {
+                if ((((i + g) >> in_qubitIndex) & 1) == result)
+                {
+                    io_psi[i + g] /= measProb;
+                }
+                else
+                {
+                    io_psi[i + g + k_range] /= measProb;
+                }    
+            }
+        }
+
+        return result;
+    }
+
+    void generateMeasureBitString(std::shared_ptr<xacc::AcceleratorBuffer> in_buffer, const std::vector<size_t>& in_bits, const KetVectorType& in_stateVec, int in_shotCount, bool in_multiThread = false)
+    {
+        if (!in_multiThread)
+        {
+            // Sequential execution
+            for (int i = 0; i < in_shotCount; ++i)
+            {
+                std::string bitString;
+                auto stateVecCopy = in_stateVec;
+                for (const auto& bit : in_bits)    
+                {
+                    bitString.append(std::to_string(applyMeasureOp(stateVecCopy, bit)));
+                }
+
+                in_buffer->appendMeasurement(bitString);
+            }
+        }
+        else
+        {
+            // Parallel execution: divide the shot counts to threads.
+            std::vector<std::string> bitStringArray(in_shotCount);
+            std::vector<std::thread> threads(getNumberOfThreads());
+            std::mutex critical;
+            for(int t = 0; t < getNumberOfThreads(); ++t)
+            {
+                threads[t] = std::thread(std::bind([&](int beginIdx, int endIdx, int threadIdx) {
+                    for(int i = beginIdx; i < endIdx; ++i)
+                    {
+                        std::string bitString;
+                        auto stateVecCopy = in_stateVec;
+                        for (const auto& bit : in_bits)    
+                        {
+                            bitString.append(std::to_string(applyMeasureOp(stateVecCopy, bit)));
+                        }
+                        bitStringArray[i] = bitString;
+                    }
+                    {
+                        // Add measurement bitstring to the buffer:
+                        std::lock_guard<std::mutex> lock(critical);
+                        for(int i = beginIdx; i < endIdx; ++i)
+                        {
+                            in_buffer->appendMeasurement(bitStringArray[i]);
+                        }
+                    }
+                }, 
+                t * in_shotCount / getNumberOfThreads(), 
+                (t+1) == getNumberOfThreads() ? in_shotCount: (t+1) * in_shotCount/getNumberOfThreads(), 
+                t));
+            }
+            std::for_each(threads.begin(),threads.end(),[](std::thread& x){
+                x.join();
+            });
+        }
+    }
+
+    // Helper to determine if shot count distribution can be simulated by 
+    // sampling the final state vector.
+    // Requirements: measure at the very end and nothing after measure.
+    bool shotCountFromFinalStateVec(const std::shared_ptr<CompositeInstruction>& in_composite)
+    {
+        InstructionIterator it(in_composite);
+        bool measureAtTheEnd = true;
+        bool measureEncountered = false;
+
+        while (it.hasNext())
+        {
+            auto nextInst = it.next();
+            if (nextInst->isEnabled())
+            {
+                if (isMeasureGate(nextInst))
+                {
+                    // Flag that we have seen a Measure gate.
+                    measureEncountered = true;
+                }
+
+                // We have seen a Measure gate but this one is not another Measure gate.
+                if (measureEncountered && !isMeasureGate(nextInst))
+                {
+                    measureAtTheEnd = false;
+                }
+            }
+        }
+
+        // If Measure gates are at the very end,
+        // this Composite can be simulated by random sampling from the state vec.
+        return measureAtTheEnd;
+    }
+}
 
 namespace xacc {
 namespace quantum {
@@ -54,12 +217,53 @@ namespace quantum {
         }
         else
         {
-           for (int i = 0; i < m_shots; ++i)
+           // Not possible to simulate shot count by direct sampling,
+           // e.g. must collapse the state vector.
+           if(!shotCountFromFinalStateVec(compositeInstruction))
            {
-                runCircuit(true);
+                for (int i = 0; i < m_shots; ++i)
+                {
+                    runCircuit(true);
+                }
+           }
+           else
+           {
+                // Index of measure bits
+                std::vector<size_t> measureBitIdxs;
+                m_visitor->initialize(buffer);
+                // Walk the IR tree, and visit each node
+                InstructionIterator it(compositeInstruction);
+                while (it.hasNext())
+                {
+                    auto nextInst = it.next();
+                    if (nextInst->isEnabled())
+                    {
+                        if (!isMeasureGate(nextInst))
+                        {
+                            nextInst->accept(m_visitor);
+                        }
+                        else
+                        {
+                            // Just collect the indices of measured qubit
+                            measureBitIdxs.emplace_back(nextInst->bits()[0]);
+                        }
+                    }
+                }
+                
+                // Run bit-string simulation
+                if (!measureBitIdxs.empty())
+                {
+                    const auto& stateVec = m_visitor->getStateVec();
+                    // Try multi-threaded execution if there are many shots.
+                    const bool multiThreadEnabled = (m_shots > 1024);
+                    generateMeasureBitString(buffer, measureBitIdxs, stateVec, m_shots, multiThreadEnabled);
+                }
+
+                m_visitor->finalize();
            }
         }
     }
+    
     void QppAccelerator::execute(std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<std::shared_ptr<CompositeInstruction>> compositeInstructions)
     {
         for (auto& f : compositeInstructions)

quantum/plugins/qpp/accelerator/QppVisitor.hpp

@@ -54,7 +54,7 @@ public:
   void visit(fSim& in_fsimGate) override;
   
   virtual std::shared_ptr<QppVisitor> clone() { return std::make_shared<QppVisitor>(); }
-
+  const KetVectorType& getStateVec() const { return m_stateVec; }
 private:
   qpp::idx xaccIdxToQppIdx(size_t in_idx) const;
   double calcExpectationValueZ() const;

quantum/plugins/qpp/tests/QppAcceleratorTester.cpp

@@ -453,6 +453,51 @@ TEST(QppAcceleratorTester, testFsim)
     }
 }
 
+TEST(QppAcceleratorTester, testDeuteronVqeH3Shots)
+{
+    // Use Qpp accelerator
+    const int nbShots = 10000;
+    auto accelerator = xacc::getAccelerator("qpp", { std::make_pair("shots", nbShots) });
+    EXPECT_EQ(accelerator->name(), "qpp");
+
+    // Create the N=3 deuteron Hamiltonian
+    auto H_N_3 = xacc::quantum::getObservable(
+        "pauli",
+        std::string("5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1 + "
+                    "9.625 - 9.625 Z2 - 3.91 X1 X2 - 3.91 Y1 Y2"));
+
+    auto optimizer = xacc::getOptimizer("nlopt");
+
+    // JIT map Quil QASM Ansatz to IR
+    xacc::qasm(R"(
+        .compiler xasm
+        .circuit deuteron_ansatz_h3_2
+        .parameters t0, t1
+        .qbit q
+        X(q[0]);
+        exp_i_theta(q, t0, {{"pauli", "X0 Y1 - Y0 X1"}});
+        exp_i_theta(q, t1, {{"pauli", "X0 Z1 Y2 - X2 Z1 Y0"}});
+    )");
+    auto ansatz = xacc::getCompiled("deuteron_ansatz_h3_2");
+
+    // Get the VQE Algorithm and initialize it
+    auto vqe = xacc::getAlgorithm("vqe");
+    vqe->initialize({std::make_pair("ansatz", ansatz),
+                    std::make_pair("observable", H_N_3),
+                    std::make_pair("accelerator", accelerator),
+                    std::make_pair("optimizer", optimizer)});
+
+    // Allocate some qubits and execute
+    auto buffer = xacc::qalloc(3);
+    vqe->execute(buffer);
+
+    // Expected result: -2.04482
+    // Tol: 0.25 (~10% of the true value)
+    // (since we are using shots, hence will introduce randomness to the optimizer)
+    std::cout << "Energy = " << (*buffer)["opt-val"].as<double>() << "\n";
+    EXPECT_NEAR((*buffer)["opt-val"].as<double>(), -2.04482, 0.25);
+}
+
 int main(int argc, char **argv) {
   xacc::Initialize();