Skip to content

Commit

Permalink
fix bugs in sample_measure
Browse files Browse the repository at this point in the history
  • Loading branch information
hhorii committed Aug 26, 2020
1 parent 0238d27 commit 986a37f
Show file tree
Hide file tree
Showing 3 changed files with 207 additions and 163 deletions.
16 changes: 0 additions & 16 deletions src/framework/utils.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -1563,22 +1563,6 @@ uint_t popcount(const uint_t count_) {
return count;
}

template <typename V, typename T>
int_t index_of(V& v, T& t, int_t first, int_t last) {
int_t mid = 0;
while (true) {
if (first >= last - 1) {
return first;
break;
}
mid = (first + last) / 2;
if (t <= v[mid])
last = mid;
else
first = mid;
}
}

//------------------------------------------------------------------------------
} // end namespace Utils
//------------------------------------------------------------------------------
Expand Down
150 changes: 3 additions & 147 deletions src/simulators/statevector/qubitvector.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,10 @@
#include <sstream>
#include <stdexcept>
#include <chrono>
#include <bitset>

#include "simulators/statevector/indexes.hpp"
#include "simulators/statevector/sampling.hpp"
#include "framework/json.hpp"
#include "framework/utils.hpp"

Expand Down Expand Up @@ -1914,152 +1916,6 @@ std::vector<double> QubitVector<data_t, Derived>::probabilities(const reg_t &qub
return probs;
}

// Get samples with saving memory with following methods of DATA.
template <typename Lambda>
reg_t sample_with_iterative_search(const std::vector<double> &rnds,
const uint_t num_qubits,
const uint_t index_size,
const Lambda probability_func,
const uint_t omp_threads) {

const int_t END = 1UL << num_qubits;
const uint_t SHOTS = rnds.size();
const int_t INDEX_SIZE = std::min(index_size, num_qubits - 1);
const int_t INDEX_END = 1UL << INDEX_SIZE;

reg_t samples;
samples.assign(SHOTS, 0);

// Initialize indices
std::vector<double> idxs;
idxs.assign(INDEX_END, 0.0);
const uint_t LOOP = (END >> INDEX_SIZE);

// Construct indices
auto idxing = [&](const int_t i)->void {
uint_t base = LOOP * i;
double total = .0;
double p = .0;
for (uint_t j = LOOP * i; j < LOOP * (i + 1); j++) {
p = probability_func(j);
total += p;
}
idxs[i] = total;
};
QV::apply_lambda(0, INDEX_END, omp_threads, idxing);

// accumulate indices
for (uint_t i = 1; i < INDEX_END; ++i)
idxs[i] += idxs[i - 1];

// reduce rounding error
double correction = 1.0 / idxs[INDEX_END - 1];
for (int_t i = 1; i < INDEX_END - 1; ++i)
idxs[i] *= correction;
idxs[INDEX_END - 1] = 1.0;

// sort random numbers
auto rnds_ = rnds;
std::sort(rnds_.begin(), rnds_.end());

// find starting index
std::vector<int_t> starts;
starts.assign(INDEX_END + 1, 0);
starts[INDEX_END] = SHOTS;

uint_t last_idx_start = 0;
for (int_t i = 1; i < INDEX_END; ++i) {
for (uint_t j = last_idx_start; j < SHOTS; ++j) {
if (rnds_[j] < idxs[i - 1])
continue;
starts[i] = j;
last_idx_start = j;
break;
}
}

// sampling
auto sampling = [&](const int_t i)->void {
uint_t start_sample_idx = starts[i];
uint_t end_sample_idx = starts[i + 1];
auto sample = LOOP * i;
double p = (i == 0) ? 0.0 : idxs[i - 1];
for (uint_t sample_idx = start_sample_idx; sample_idx < end_sample_idx; ++sample_idx) {
auto rnd = rnds_[sample_idx];
for (; sample < LOOP * (i + 1) - 1; ++sample) {
p += probability_func(sample);
if (rnd < p){
break;
}
}
samples[sample_idx] = sample;
}
};
QV::apply_lambda(0, INDEX_END, omp_threads, sampling);

return samples;
}

// Get samples with consumption of memory with following methods of DATA.
template <typename Lambda>
reg_t sample_with_binary_search(const std::vector<double> &rnds,
const uint_t num_qubits,
const Lambda probability_func,
const uint_t omp_threads) {

const int_t END = 1UL << num_qubits;
const uint_t SHOTS = rnds.size();
const uint_t PARTITION_SIZE = num_qubits < 15 ? 0: 10;
const int_t PARTITION_END = 1UL << PARTITION_SIZE;
const uint_t LOOP = (END >> PARTITION_SIZE);

reg_t samples;
samples.assign(SHOTS, 0);

std::vector<std::vector<double>> acc_probs_list;
std::vector<reg_t> acc_idxs_list;
std::vector<double> start_probs;
std::vector<double> end_probs;

acc_probs_list.assign(PARTITION_END, std::vector<double>());
acc_idxs_list.assign(PARTITION_END, reg_t());
start_probs.assign(PARTITION_END, 0.);
end_probs.assign(PARTITION_END, 0.);

// inclusive scan for each partition
auto scan = [&](const int_t i)->void {
double accumulated = .0;
for (uint_t j = LOOP * i; j < LOOP * (i + 1); j++) {
auto norm = probability_func(j);
if (!AER::Linalg::almost_equal(norm, 0.0)) {
accumulated += norm;
acc_probs_list[i].push_back(accumulated);
acc_idxs_list[i].push_back(j);
}
}
end_probs[i] = accumulated;
};
QV::apply_lambda(0, PARTITION_END, omp_threads, scan);

for (int_t i = 1; i < PARTITION_END; ++i)
start_probs[i] = end_probs[i -1] + start_probs[i - 1];

// sampling
auto sampling = [&](const int_t i)->void {
double rnd = rnds[i];

// binary search for partition
int_t partition_idx = AER::Utils::index_of(start_probs, rnd, 0, PARTITION_END - 1);
rnd -= start_probs[partition_idx];
// binary search for which range rnd is in
int_t sample_idx = AER::Utils::index_of(acc_probs_list[partition_idx], rnd, 0, acc_probs_list[partition_idx].size() - 1);
samples[i] = acc_idxs_list[partition_idx][sample_idx];
};
QV::apply_lambda(0, SHOTS, omp_threads, sampling);

return samples;
}

template <typename data_t, typename Derived>
reg_t QubitVector<data_t, Derived>::sample_measure(const std::vector<double> &rnds) const {

Expand All @@ -2072,8 +1928,8 @@ reg_t QubitVector<data_t, Derived>::sample_measure(const std::vector<double> &rn
ret = sample_with_binary_search(rnds, num_qubits(), lambda, omp_threads_managed());
else
ret = sample_with_iterative_search(rnds, num_qubits(), sample_measure_index_size_, lambda, omp_threads_managed());
return ret;

return ret;
}

/*******************************************************************************
Expand Down
Loading

0 comments on commit 986a37f

Please sign in to comment.