Skip to content

Commit

Permalink
Fix measure in stabilizer method (#1895)
Browse files Browse the repository at this point in the history 
* Fix measure in stabilizer

* add test case for measuring stabilizer

---------

Co-authored-by: Hiroshi Horii <hhorii@users.noreply.github.com>
  • Loading branch information
@doichanj @hhorii
doichanj and hhorii authored Sep 7, 2023
1 parent 6460e49 commit 4fb99dd
Show file tree
Hide file tree
Showing 3 changed files with 137 additions and 71 deletions.
4 changes: 4 additions & 0 deletions releasenotes/notes/fix_stabilizer_measure-49f9b4261dfaa4d3.yaml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
---
fixes:
- |
This release fixes an issue in measurement function of stabilizer simulator
141 changes: 70 additions & 71 deletions src/simulators/stabilizer/clifford.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -422,7 +422,6 @@ bool Clifford::measure_and_update(const uint64_t qubit,
auto anticom = z_anticommuting(qubit);

int nid = omp_get_num_threads();

if (anticom.first) {
bool outcome = (randint == 1);
auto row = anticom.second;
Expand Down Expand Up @@ -515,90 +514,86 @@ bool Clifford::measure_and_update(const uint64_t qubit,
return outcome;
} else {
// Deterministic outcome
bool outcome = false;
uint_t outcome = 0;
Pauli::Pauli<BV::BinaryVector> accum(num_qubits_);
uint64_t blocks = destabilizer_phases_.blockLength();
uint_t blocks = destabilizer_phases_.blockLength();

if (blocks < 2) {
for (uint64_t i = 0; i < num_qubits_; i++) {
if (destabilizer_table_[qubit].X[i]) {
bool b0 = false, b1 = false;
for (size_t q = 0; q < num_qubits_; q++) {
bool t0, t1, add;
bool accumX = accum.X[q];
bool accumZ = accum.Z[q];

t0 = accumX & stabilizer_table_[q].Z[i];
t1 = accumZ ^ stabilizer_table_[q].X[i];

b1 ^= (t0 & b0);
b0 ^= t0;
b1 ^= (t0 & t1);

t0 = stabilizer_table_[q].X[i] & accumZ;
t1 = stabilizer_table_[q].Z[i] ^ accumX;
t1 ^= t0;

b1 ^= (t0 & b0);
b0 ^= t0;
b1 ^= (t0 & t1);

accum.X.setValue(stabilizer_table_[q].X[i] ^ accum.X[q], q);
accum.Z.setValue(stabilizer_table_[q].Z[i] ^ accum.Z[q], q);
}
b1 ^= (stabilizer_phases_[i] ^ outcome);

if (b0) {
throw std::runtime_error("Clifford: rowsum error");
}
outcome = b1;
for (uint_t ib = 0; ib < blocks; ib++) {
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 {
uint64_t blockSize = destabilizer_phases_.blockSize();
uint_t blockSize = destabilizer_phases_.blockSize();

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

uint64_t exponent_l = 0;
uint64_t exponent_lc = 0;
uint64_t exponent_h = 0;
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) {
uint64_t qs = qq * blockSize;
uint64_t qe = qs + blockSize;
uint_t qs = qq * blockSize;
uint_t qe = qs + blockSize;
if (qe > num_qubits_)
qe = num_qubits_;

uint64_t local_exponent_l = 0;
uint64_t local_exponent_h = 0;

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

// set accum for this block
uint64_t accumX = destabilizer_mask & sX;
uint64_t accumZ = destabilizer_mask & sZ;
for (int b = 1; b < blockSize; b *= 2) {
accumX ^= (accumX << b);
accumZ ^= (accumZ << b);
}
accumX ^= (0ull - (uint64_t)accum.X[q]);
accumZ ^= (0ull - (uint64_t)accum.Z[q]);
accum.X.setValue((accumX >> (blockSize - 1)), q);
accum.Z.setValue((accumZ >> (blockSize - 1)), q);
for (uint_t q = qs; q < qe; q++) {
uint_t sX = stabilizer_table_[q].X(ii);
uint_t sZ = stabilizer_table_[q].Z(ii);

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

// exponents for this block
uint64_t t0, t1;
uint_t t0, t1;

t0 = accumX & sZ;
t1 = accumZ ^ sX;
Expand All @@ -614,6 +609,12 @@ bool Clifford::measure_and_update(const uint64_t qubit,
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
Expand All @@ -627,16 +628,14 @@ bool Clifford::measure_and_update(const uint64_t qubit,
(num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0,
blocks, measure_determinisitic_func, omp_threads_);

exponent_h ^=
(exponent_lc ^
exponent_l); // if exponent_l is 0 and any of local_exponent_l is
// 1, then flip exponent_h

exponent_h ^= (stabilizer_phases_(ii) & destabilizer_mask);
outcome ^= ((AER::Utils::popcount(exponent_h) & 1) != 0);
// 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);
}
}
return outcome;
return ((AER::Utils::popcount(outcome) & 1) != 0);
}
}

Expand Down
63 changes: 63 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 @@ -23,6 +23,7 @@
from qiskit.circuit.library import QuantumVolume
from qiskit.quantum_info.random import random_unitary
from test.terra.backends.simulator_test_case import SimulatorTestCase, supported_methods
import numpy as np

SUPPORTED_METHODS = [
"automatic",
Expand Down Expand Up @@ -199,6 +200,68 @@ def test_measure_nondeterministic_multi_qubit_without_sampling(self, method, dev
self.compare_counts(result, circuits, targets, delta=delta * shots)
self.compare_result_metadata(result, circuits, "measure_sampling", False)

# ---------------------------------------------------------------------
# Test stabilizer measure
# ---------------------------------------------------------------------
@supported_methods(["stabilizer"])
def test_measure_stablizer_64bit(self, method, device):
backend = self.backend(method=method, device=device)
shots = 10000
delta = 0.05
circ = QuantumCircuit(65, 32)

circ.reset(0)
for i in range(0, 30, 6):
circ.h(i)
circ.h(i + 4)
circ.h(30)
circ.h(31)

for i in range(1, 32, 2):
circ.cx(i + 32, i)
for i in range(0, 30, 6):
circ.cx(i, i + 32)
circ.cx(i + 4, i + 36)
circ.cx(30, 62)

for i in range(1, 30, 2):
circ.cx(i + 35, i)
for i in range(4, 32, 4):
circ.cx(i, i + 29)

for i in range(0, 30, 2):
circ.cx(i + 35, i)
for i in range(1, 30, 6):
circ.cx(i, i + 33)
circ.cx(i + 2, i + 35)
circ.cx(31, 64)

for i in range(0, 32):
circ.measure(i, i)
result = backend.run(circ, shots=shots).result()
counts = result.get_counts()
self.assertSuccess(result)

n_anc = 32
totals = np.zeros(n_anc, dtype=int)
for outcomes, num_counts in counts.items():
new_totals = num_counts * np.array([int(bit) for bit in outcomes][::-1])
assert len(new_totals) == n_anc
totals += new_totals
output = {}
for i in range(0, 32):
output[hex(i)] = totals[i]

targets = {}
for i in range(0, 30, 3):
targets[hex(i)] = shots / 2
targets[hex(i + 1)] = shots / 2
targets[hex(i + 2)] = 0
targets[hex(30)] = shots / 2
targets[hex(31)] = shots / 2

self.assertDictAlmostEqual(output, targets, delta=delta * shots)

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

0 comments on commit 4fb99dd

Please sign in to comment.