Skip to content

Commit

Permalink
Fix deterministic measure of stabilizer (Qiskit#2132)
Browse files Browse the repository at this point in the history 
* Fix deterministic measure of stabilizer

* fix accumulation in loop

* format

---------

Co-authored-by: Hiroshi Horii <hhorii@users.noreply.github.com>
  • Loading branch information
@doichanj @hhorii
doichanj and hhorii committed May 28, 2024
1 parent 6281a23 commit 42c82c1
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 118 deletions.
5 changes: 5 additions & 0 deletions releasenotes/notes/fix_stabilizer_measure-a06d761eba2546d2.yaml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
---
fixes:
- |
Fix of deterministic measure of stabilizer and change of OpenMP
parallelization of loop
186 changes: 68 additions & 118 deletions src/simulators/stabilizer/clifford.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -456,9 +456,8 @@ bool Clifford::measure_and_update(const uint64_t qubit,

t0 = rX & destabilizer_table_[q].Z(i);
t1 = rZ ^ destabilizer_table_[q].X(i);
t1 ^= t0;

d1 ^= (t0 & d0);
d1 ^= (t0 & (~d0));
d0 ^= t0;
d1 ^= (t0 & t1);

Expand All @@ -475,9 +474,8 @@ bool Clifford::measure_and_update(const uint64_t qubit,

t0 = rX & stabilizer_table_[q].Z(i);
t1 = rZ ^ stabilizer_table_[q].X(i);
t1 ^= t0;

s1 ^= (t0 & s0);
s1 ^= (t0 & (~s0));
s0 ^= t0;
s1 ^= (t0 & t1);

Expand Down Expand Up @@ -515,127 +513,79 @@ bool Clifford::measure_and_update(const uint64_t qubit,
} else {
// Deterministic outcome
uint_t outcome = 0;
Pauli::Pauli<BV::BinaryVector> accum(num_qubits_);
uint_t blocks = destabilizer_phases_.blockLength();

if (blocks < 2) {
auto measure_determinisitic_func = [this, blocks, qubit](AER::int_t q) {
uint_t accumX = 0;
uint_t accumZ = 0;
uint_t exponent_l = 0ull;
uint_t exponent_h = 0ull;
bool accumX_prev = false;
bool accumZ_prev = false;

for (uint_t ib = 0; ib < blocks; ib++) {
uint_t tl, th, add;
uint_t destabilizer_mask = destabilizer_table_[qubit].X(ib);
uint_t exponent_l = 0ull;
uint_t exponent_h = 0ull;

for (uint_t q = 0; q < num_qubits_; q++) {
uint_t tl, th, add;
uint_t accumX = 0ull - (uint_t)accum.X[q];
uint_t accumZ = 0ull - (uint_t)accum.Z[q];

tl = accumX & stabilizer_table_[q].Z(ib);
th = accumZ ^ stabilizer_table_[q].X(ib);

add = tl & exponent_l;
exponent_l ^= tl;
exponent_h ^= add;
exponent_h ^= (tl & th);

tl = stabilizer_table_[q].X(ib) & accumZ;
th = stabilizer_table_[q].Z(ib) ^ accumX;
th ^= tl;

add = tl & exponent_l;
exponent_l ^= tl;
exponent_h ^= add;
exponent_h ^= (tl & th);

add = stabilizer_table_[q].X(ib) & destabilizer_mask;
accumX &= AER::Utils::popcount(add) & 1;
add = stabilizer_table_[q].Z(ib) & destabilizer_mask;
accumZ &= AER::Utils::popcount(add) & 1;

accum.X.setValue((bool)accumX, q);
accum.Z.setValue((bool)accumZ, q);
}
exponent_h ^= stabilizer_phases_(ib);
outcome ^= (exponent_h & destabilizer_mask);

if ((exponent_l & destabilizer_mask) != 0) {
throw std::runtime_error("Clifford: rowsum error");
}
}
} else {
uint_t blockSize = destabilizer_phases_.blockSize();

// loop for cache blocking
for (uint_t ii = 0; ii < blocks; ii++) {
uint_t destabilizer_mask = destabilizer_table_[qubit].X(ii);
if (destabilizer_mask == 0)
continue;

uint_t exponent_l = 0;
uint_t exponent_lc = 0;
uint_t exponent_h = 0;

auto measure_determinisitic_func =
[this, &accum, &exponent_l, &exponent_lc, &exponent_h, blocks,
blockSize, destabilizer_mask, ii](AER::int_t qq) {
uint_t qs = qq * blockSize;
uint_t qe = qs + blockSize;
if (qe > num_qubits_)
qe = num_qubits_;

uint_t local_exponent_l = 0;
uint_t local_exponent_h = 0;

for (uint_t q = qs; q < qe; q++) {
uint_t sX = stabilizer_table_[q].X(ii);
uint_t sZ = stabilizer_table_[q].Z(ii);

uint_t accumX = (0ull - (uint_t)accum.X[q]);
uint_t accumZ = (0ull - (uint_t)accum.Z[q]);

// exponents for this block
uint_t t0, t1;

t0 = accumX & sZ;
t1 = accumZ ^ sX;

local_exponent_h ^= (t0 & local_exponent_l);
local_exponent_l ^= t0;
local_exponent_h ^= (t0 & t1);

t0 = sX & accumZ;
t1 = sZ ^ accumX;
t1 ^= t0;

local_exponent_h ^= (t0 & local_exponent_l);
local_exponent_l ^= t0;
local_exponent_h ^= (t0 & t1);

// update accum
accumX &= AER::Utils::popcount(sX & destabilizer_mask) & 1;
accum.X.setValue((accumX != 0), q);
accumZ &= AER::Utils::popcount(sZ & destabilizer_mask) & 1;
accum.Z.setValue((accumZ != 0), q);
}

#pragma omp atomic
exponent_lc |= local_exponent_l;
#pragma omp atomic
exponent_l ^= local_exponent_l;
#pragma omp atomic
exponent_h ^= local_exponent_h;
};
AER::Utils::apply_omp_parallel_for(
(num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0,
blocks, measure_determinisitic_func, omp_threads_);

// if exponent_l is 0 and any of local_exponent_l is
// 1, then flip exponent_h
exponent_h ^= (exponent_lc ^ exponent_l);
exponent_h ^= stabilizer_phases_(ii);
outcome ^= (exponent_h & destabilizer_mask);
uint_t sX = stabilizer_table_[q].X(ib) & destabilizer_mask;
uint_t sZ = stabilizer_table_[q].Z(ib) & destabilizer_mask;

// accumulate for 64 bits block
accumX = sX ^ (uint_t)accumX_prev;
accumZ = sZ ^ (uint_t)accumZ_prev;
accumX ^= (accumX << 1);
accumZ ^= (accumZ << 1);
accumX ^= (accumX << 2);
accumZ ^= (accumZ << 2);
accumX ^= (accumX << 4);
accumZ ^= (accumZ << 4);
accumX ^= (accumX << 8);
accumZ ^= (accumZ << 8);
accumX ^= (accumX << 16);
accumZ ^= (accumZ << 16);
accumX ^= (accumX << 32);
accumZ ^= (accumZ << 32);
// store for next iteration
accumX_prev = ((accumX >> 63) & 1) != 0;
accumZ_prev = ((accumZ >> 63) & 1) != 0;
// correct for this iteration
accumX ^= sX;
accumZ ^= sZ;
accumX &= destabilizer_mask;
accumZ &= destabilizer_mask;

tl = accumX & sZ;
th = accumZ ^ sX;

add = tl & exponent_l;
exponent_l ^= tl;
exponent_h ^= add;
exponent_h ^= (tl & th);

tl = sX & accumZ;
th = sZ ^ accumX;

add = tl & (~exponent_l);
exponent_l ^= tl;
exponent_h ^= add;
exponent_h ^= (tl & th);
}
// convert 2-bits x 64 integer into bit count here
return AER::Utils::popcount(exponent_h) * 2 +
AER::Utils::popcount(exponent_l);
};

outcome = AER::Utils::apply_omp_parallel_for_reduction_int(
(num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0,
num_qubits_, measure_determinisitic_func, omp_threads_);

uint_t stab_h = 0ull;
for (uint_t ib = 0; ib < blocks; ib++) {
stab_h ^= (destabilizer_table_[qubit].X(ib) & stabilizer_phases_(ib));
}
return ((AER::Utils::popcount(outcome) & 1) != 0);
outcome += AER::Utils::popcount(stab_h) * 2;

return ((outcome & 3) == 2);
}
}

Expand Down
32 changes: 32 additions & 0 deletions test/terra/backends/aer_simulator/test_measure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,38 @@ def test_measure_sampling_large_ghz_stabilizer(self, method, device, num_qubits)
targets["1" * num_qubits] = shots / 2
self.assertDictAlmostEqual(counts, targets, delta=delta * shots)

@supported_methods(["stabilizer"])
def test_measure_stablizer_deterministic(self, method, device):
"""Test stabilizer measure for deterministic case"""
backend = self.backend(method=method, device=device)
shots = 10000
qc = QuantumCircuit(5, 1)
qc.h([2, 4])
qc.cx(2, 0)
qc.s(0)
qc.cx(4, 2)
qc.h(0)
qc.cx(2, 3)
qc.s(4)
qc.cx(1, 0)
qc.h([3, 4])
qc.cx(3, 2)
qc.h(3)
qc.cx(0, 3)
qc.cx(3, 1)
qc.s(0)
qc.s(1)
qc.h(0)
qc.s(0)
qc.cx(4, 0)
qc.cx(0, 1)
qc.measure(1, 0)
result = backend.run(qc, shots=shots).result()
counts = result.get_counts()
self.assertSuccess(result)

self.assertDictEqual(counts, {"1": shots})

# ---------------------------------------------------------------------
# Test MPS algorithms for measure
# ---------------------------------------------------------------------
Expand Down

0 comments on commit 42c82c1

Please sign in to comment.