In Stim, after performing syndrome measurements circuit using stabilizer codes like surface codes, how can we know the expected value of errors occurring on each data qubit? Even if there are no errors in the circuit, measuring each data qubit transversally after the syndrome measurement will result in the measurement outcome flipping non-deterministically (up to the stabilizer). This means that simply sampling the results of the transversal measurement of each data qubit to obtain the expected value won't yield the correct expectation value. I'd like to know if there's another way.
1 Answer
$\begingroup$
$\endgroup$
Since stim v1.12.0, you can use stim.FlipSimulator to determine the rate at which data qubits are flipped (or whether, in a specific shot, it has been flipped a particular way).
Truncate the circuit to just before the data measurements, run the circuit in the flip simulator (with stabilizer randomization disabled), and count how often you see each Pauli.
import numpy as np
import stim
# Get a surface code circuit.
circuit = stim.Circuit.generated(
"surface_code:rotated_memory_x",
distance=5,
rounds=15,
after_clifford_depolarization=1e-3,
before_round_data_depolarization=1e-3,
before_measure_flip_probability=1e-3,
after_reset_flip_probability=1e-3,
)
# Truncate the circuit so it stops just before the data measurements.
last_measurement_layer = len(circuit) - 1
while circuit[last_measurement_layer].name != 'MR':
last_measurement_layer -= 1
circuit = circuit[:last_measurement_layer]
# Collect stats over a million+ shots, in batches of 1024 (takes ~10 seconds).
i_hits = np.zeros(shape=circuit.num_qubits, dtype=np.uint64)
x_hits = np.zeros(shape=circuit.num_qubits, dtype=np.uint64)
y_hits = np.zeros(shape=circuit.num_qubits, dtype=np.uint64)
z_hits = np.zeros(shape=circuit.num_qubits, dtype=np.uint64)
for _ in range(1000):
sim = stim.FlipSimulator(
batch_size=1024,
disable_stabilizer_randomization=True,
)
sim.do(circuit)
# Count number of times each Pauli occurred on each qubit.
instance_paulis: stim.PauliString
for instance_paulis in sim.peek_pauli_flips():
xs, zs = instance_paulis.to_numpy(bit_packed=False)
i_hits += ~xs & ~zs
x_hits += xs & ~zs
y_hits += xs & zs
z_hits += ~xs & zs
# Print results.
qubit_coords = circuit.get_final_qubit_coordinates()
for q, coords in qubit_coords.items():
i, x, y, z = i_hits[q], x_hits[q], y_hits[q], z_hits[q]
t = i + x + y + z
print(f"qubit at {str(coords):>14}: X={x/t:<16} Y={y/t:<16} Z={z/t:<16}")
Example results (note this includes data and measure qubits):
qubit at [1.0, 1.0]: X=0.00206640625 Y=0.0008642578125 Z=0.0033662109375
qubit at [2.0, 0.0]: X=0.006416015625 Y=0.0011796875 Z=0.0012509765625
qubit at [3.0, 1.0]: X=0.0041142578125 Y=0.0010810546875 Z=0.0027060546875
qubit at [5.0, 1.0]: X=0.004029296875 Y=0.0011611328125 Z=0.003640625
qubit at [6.0, 0.0]: X=0.006451171875 Y=0.001201171875 Z=0.0011552734375
qubit at [7.0, 1.0]: X=0.0041064453125 Y=0.001087890625 Z=0.00265625
qubit at [9.0, 1.0]: X=0.00259765625 Y=0.00084765625 Z=0.0024365234375
qubit at [1.0, 3.0]: X=0.0034638671875 Y=0.0011806640625 Z=0.003703125
qubit at [2.0, 2.0]: X=0.007291015625 Y=0.00102734375 Z=0.0010830078125
qubit at [3.0, 3.0]: X=0.004375 Y=0.0014482421875 Z=0.0039501953125
qubit at [4.0, 2.0]: X=0.0123115234375 Y=0.0018193359375 Z=0.0017705078125
qubit at [5.0, 3.0]: X=0.0043759765625 Y=0.0014013671875 Z=0.0039853515625
qubit at [6.0, 2.0]: X=0.0089931640625 Y=0.001046875 Z=0.0010712890625
qubit at [7.0, 3.0]: X=0.004392578125 Y=0.0014306640625 Z=0.0039423828125
qubit at [8.0, 2.0]: X=0.0118388671875 Y=0.00170703125 Z=0.0017861328125
qubit at [9.0, 3.0]: X=0.0018125 Y=0.0011689453125 Z=0.002677734375
qubit at [10.0, 2.0]: X=0.0059580078125 Y=0.000498046875 Z=0.00056640625
qubit at [0.0, 4.0]: X=0.003765625 Y=0.000537109375 Z=0.0005283203125
qubit at [1.0, 5.0]: X=0.002279296875 Y=0.0011962890625 Z=0.0036416015625
qubit at [2.0, 4.0]: X=0.0132919921875 Y=0.0018017578125 Z=0.0017685546875
qubit at [3.0, 5.0]: X=0.004310546875 Y=0.001388671875 Z=0.00394140625
qubit at [4.0, 4.0]: X=0.0094189453125 Y=0.00109765625 Z=0.0010927734375
qubit at [5.0, 5.0]: X=0.004279296875 Y=0.001435546875 Z=0.0039130859375
qubit at [6.0, 4.0]: X=0.01394921875 Y=0.0017841796875 Z=0.001677734375
qubit at [7.0, 5.0]: X=0.0043056640625 Y=0.001427734375 Z=0.0039638671875
qubit at [8.0, 4.0]: X=0.0096455078125 Y=0.0010634765625 Z=0.0010810546875
qubit at [9.0, 5.0]: X=0.0028818359375 Y=0.0011396484375 Z=0.002625
qubit at [1.0, 7.0]: X=0.003470703125 Y=0.0011318359375 Z=0.0036708984375
qubit at [2.0, 6.0]: X=0.00767578125 Y=0.0010361328125 Z=0.0010830078125
qubit at [3.0, 7.0]: X=0.0042568359375 Y=0.00137109375 Z=0.0040126953125
qubit at [4.0, 6.0]: X=0.0139052734375 Y=0.0017578125 Z=0.0016982421875
qubit at [5.0, 7.0]: X=0.004298828125 Y=0.001373046875 Z=0.0039921875
qubit at [6.0, 6.0]: X=0.009423828125 Y=0.001046875 Z=0.00103515625
qubit at [7.0, 7.0]: X=0.0043837890625 Y=0.001388671875 Z=0.003896484375
qubit at [8.0, 6.0]: X=0.013353515625 Y=0.001693359375 Z=0.0017841796875
qubit at [9.0, 7.0]: X=0.001828125 Y=0.00116796875 Z=0.0025703125
qubit at [10.0, 6.0]: X=0.0066552734375 Y=0.0005341796875 Z=0.000509765625
qubit at [0.0, 8.0]: X=0.003779296875 Y=0.0005263671875 Z=0.00051953125
qubit at [1.0, 9.0]: X=0.0020869140625 Y=0.00085546875 Z=0.001837890625
qubit at [2.0, 8.0]: X=0.0133505859375 Y=0.0017529296875 Z=0.0017900390625
qubit at [3.0, 9.0]: X=0.0029619140625 Y=0.0011484375 Z=0.00323828125
qubit at [4.0, 8.0]: X=0.008765625 Y=0.001076171875 Z=0.0010576171875
qubit at [5.0, 9.0]: X=0.0030146484375 Y=0.001171875 Z=0.0021728515625
qubit at [6.0, 8.0]: X=0.0140751953125 Y=0.0018154296875 Z=0.00168359375
qubit at [7.0, 9.0]: X=0.0030009765625 Y=0.001162109375 Z=0.0032080078125
qubit at [8.0, 8.0]: X=0.0088955078125 Y=0.00108984375 Z=0.0010390625
qubit at [9.0, 9.0]: X=0.0015361328125 Y=0.0008759765625 Z=0.001861328125
qubit at [4.0, 10.0]: X=0.009263671875 Y=0.0011455078125 Z=0.0011689453125
qubit at [8.0, 10.0]: X=0.00867578125 Y=0.0012060546875 Z=0.001201171875
