Commit 39d4801d authored by Thien Nguyen's avatar Thien Nguyen
Browse files

Rework bitstring sampling from reconstructed state-vector (small number of qubits)

Batch samping whole bitstrings (all shots) at once.
parent df4a8ca3
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -113,7 +113,7 @@ public:

    if (config.keyExists<int>("seed")) {
      const auto seed = config.get<int>("seed");
      randomEngine::setSeed(seed);
      randomEngine::get_instance().setSeed(seed);
    }

    // Updated the cached configurations (to be sent on to visitor)
+33 −0
Original line number Diff line number Diff line
@@ -122,6 +122,39 @@ void ApplyCNOTGate(StateVectorType& io_psi, size_t in_controlIndex, size_t in_ta
  }
}

template <typename ElementType, typename IndexType>
std::vector<std::string>
GenerateSamples(const std::vector<ElementType> &state, uint64_t num_samples,
                const std::vector<IndexType> &meas_bits) {
  const auto toBitString = [&](uint64_t val) -> std::string {
    std::string result(meas_bits.size(), '0');
    for (size_t i = 0; i < meas_bits.size(); ++i) {
      if (val & (1ULL << meas_bits[i])) {
        result[i] = '1';
      }
    }
    return result;
  };
  std::vector<std::string> bitstrings;
  if (num_samples > 0) {
    const uint64_t size = state.size();
    const auto rs =
        tnqvm::randomEngine::get_instance().sortedRandProbs(num_samples);
    bitstrings.reserve(num_samples);
    double csum = 0.0;
    uint64_t m = 0;
    for (uint64_t k = 0; k < size; ++k) {
      csum += std::norm(state[k]);
      while (rs[m] < csum && m < num_samples) {
        bitstrings.emplace_back(toBitString(k));
        ++m;
      }
    }
  }

  return bitstrings;
}

template<typename ElementType>
bool ApplyMeasureOp(std::vector<ElementType>& io_psi, size_t in_qubitIndex)
{
+31 −16
Original line number Diff line number Diff line
#pragma once
#include <algorithm>
#include <random>
#include <mutex>
#include <utility>
#include <vector>
namespace tnqvm {
struct randomEngine {
  randomEngine(const randomEngine &) = delete;
  randomEngine &operator=(const randomEngine &) = delete;
  double randProb() {
    const auto val = std::uniform_real_distribution<double>(0.0, 1.0)(m_engine);
    return val;
    std::lock_guard<std::mutex> lock(m_mutex);
    return std::uniform_real_distribution<double>(0.0, 1.0)(m_engine);
  }

  thread_local static randomEngine &get_instance() {
    thread_local static randomEngine instance;
    thread_local static size_t seed = 0;
    if (seed != globalSeed) {
      instance.m_engine.seed(globalSeed);
      seed = globalSeed;
    }
  static randomEngine &get_instance() {
    static randomEngine instance;
    return instance;
  }
  std::mt19937_64 m_engine;
  static inline size_t globalSeed = []() {
    std::random_device rd;
    return rd();
  }();
  static void setSeed(size_t seed) { globalSeed = seed; }

  void setSeed(size_t seed) {
    std::lock_guard<std::mutex> lock(m_mutex);
    m_engine.seed(seed);
  }

  std::vector<double> sortedRandProbs(uint64_t num_samples) {
    std::vector<double> rs;
    rs.reserve(num_samples + 1);
    std::lock_guard<std::mutex> lock(m_mutex);
    for (uint64_t i = 0; i < num_samples; ++i) {
      rs.emplace_back(
          std::uniform_real_distribution<double>(0.0, 1.0)(m_engine));
    }
    std::sort(rs.begin(), rs.end());
    return rs;
  }

private:
  randomEngine() = default;
  randomEngine() {
    std::random_device rd;
    setSeed(rd());
  }
  std::mt19937_64 m_engine;
  std::mutex m_mutex;
};
} // namespace tnqvm
 No newline at end of file
+6 −56
Original line number Diff line number Diff line
@@ -259,7 +259,7 @@ void ExatnMpsVisitor::initialize(std::shared_ptr<AcceleratorBuffer> buffer, int
    if (options.keyExists<double>("svd-cutoff"))
    {
        m_svdCutoff = options.get<double>("svd-cutoff");
        std::cout << "[DEBUG] SVD Cut-off = " << m_svdCutoff << "\n";
        // std::cout << "[DEBUG] SVD Cut-off = " << m_svdCutoff << "\n";
    }

    // Max bond dimension: take precedent over svd-cutoff
@@ -267,7 +267,7 @@ void ExatnMpsVisitor::initialize(std::shared_ptr<AcceleratorBuffer> buffer, int
    if (options.keyExists<int>("max-bond-dim"))
    {
        m_maxBondDim = options.get<int>("max-bond-dim");
        std::cout << "[DEBUG] Max bond dimension = " << m_maxBondDim << "\n";
        // std::cout << "[DEBUG] Max bond dimension = " << m_maxBondDim << "\n";
    }

    m_buffer = std::move(buffer);
@@ -2202,59 +2202,9 @@ void ExatnMpsVisitor::evaluateTensorNetwork(exatn::numerics::TensorNetwork& io_t

void ExatnMpsVisitor::addMeasureBitStringProbability(const std::vector<size_t>& in_bits, const std::vector<std::complex<double>>& in_stateVec, int in_shotCount)
{
    // Factor to determine if we should spawn threads to simulate bitstring sampling.
    const int MULT_FACTOR = 100;
    if (getNumberOfThreads() < 2 || in_shotCount < MULT_FACTOR * getNumberOfThreads())
    {
        // 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)));
            }

            m_buffer->appendMeasurement(bitString);
        }
    }
    else
    {
        // Parallel execution
        std::vector<std::string> bitStringArray(in_shotCount);
        assert(bitStringArray.size() == 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)
                    {
                        m_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();
        });
  const auto measBitStrs = GenerateSamples(in_stateVec, in_shotCount, in_bits);
  for (const auto &sample : measBitStrs) {
    m_buffer->appendMeasurement(sample);
  }
}

+1 −1
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ TEST(MpsMeasurementTester, checkSimple)
TEST(MpsMeasurementTester, checkRandomSeed) 
{
  std::vector<std::map<std::string, int>> results;
  constexpr int NB_TESTS = 4;
  constexpr int NB_TESTS = 20;
  for (int i = 0; i < NB_TESTS; i++) {
    auto xasmCompiler = xacc::getCompiler("xasm");
    auto ir = xasmCompiler->compile(R"(__qpu__ void test1(qbit q) {
Loading