<jupyter_start><jupyter_text>Copyright Rigetti Computing 2021.  Sampling and Readout Errors using pyQuil Import relevant modules.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

import itertools

import numpy as np

from pyquil import get_qc, Program
from pyquil.api import WavefunctionSimulator
from pyquil.gates import *
from pyquil.paulis import *

from tqdm import tqdm

import networkx as nx<jupyter_output><empty_output><jupyter_text>Sampling Error Previously, we used the `WavefunctionSimulator`. This time, we will primarily use a `QuantumComputer`.<jupyter_code>wfn_sim = WavefunctionSimulator()

n_qubits = 2
qc = get_qc(f'{n_qubits}q-qvm')
qubits = qc.qubits()
print(qubits)<jupyter_output>[0, 1]<jupyter_text>Let us first compute what the expectation value of the observable $Z$ is in a fair (quantum) dice: the Hadamard state.<jupyter_code>prog = Program(H(0))
Z_expectation = wfn_sim.expectation(prog, sZ(0)+0.0)
print(Z_expectation)<jupyter_output>0j<jupyter_text>Now let us try to estimate it via sampling. First, we simply obtain some shots.<jupyter_code>n_shots = 1000

prog = Program()
ro = prog.declare('ro', 'BIT', 1)
prog += H(0)
prog += MEASURE(0, ro[0])
prog.wrap_in_numshots_loop(n_shots)

exe = qc.compile(prog)

samples = qc.run(exe)
print(samples.shape)
print(samples[:5])<jupyter_output>(1000, 1)
[[1]
 [1]
 [0]
 [1]
 [0]]<jupyter_text>Next, we estimate the expectation value $\langle Z \rangle$. Remember that $\langle 0 \vert Z \vert 0 \rangle=1$ and $\langle 1 \vert Z \vert 1 \rangle=-1$.<jupyter_code>Z_eigs = 1 - 2*samples
Z_expect = np.mean(Z_eigs)
print(Z_expect)<jupyter_output>-0.01<jupyter_text>Ok, so it's close, though not exactly equal, to our ideal value. Let us perform several of the above experiments, and obtain a distribution of such averages. (Recall: The Central Limit Theorem says this should look like a bell-shaped curve).<jupyter_code>n_expts = 1000

Z_expects = []

for _ in tqdm(range(n_expts)):
    
    samples = qc.run(exe)
    Z_eigs = 1 - 2*samples
    Z_expect = np.mean(Z_eigs)
    Z_expects.append(Z_expect)
    
plt.figure(figsize=(12, 8))
plt.hist(Z_expects, bins=50)
plt.show()<jupyter_output>100%|██████████| 1000/1000 [00:16<00:00, 59.35it/s]<jupyter_text>We can also estimate expectation values of other observables, such as $X$ or $Y$. But in order to do so, we need to rotate into the appropriate basis, since our quantum computer can really only measure in the $Z$ (computational) basis.In order to measure $X$, note that $HZH = X$, so that $\langle \psi \vert X \vert \psi \rangle = \langle \psi \vert HZH \vert \psi \rangle$, so that we can simply apply an $H$ gate before measuring, then carrying out the usual procedure. First, let's compute what the ideal value is.<jupyter_code>prog = Program(H(0))
X_expectation = wfn_sim.expectation(prog, sX(0)+0.0)
print(X_expectation)<jupyter_output>(0.9999999999999998+0j)<jupyter_text>Next, let's estimate via sampling.<jupyter_code>n_shots = 1000

prog = Program()
ro = prog.declare('ro', 'BIT', 1)
prog += H(0)

# apply (another) H to perform a measurement in the X basis
prog += H(0)

# now measure as usual
prog += MEASURE(0, ro[0])
prog.wrap_in_numshots_loop(n_shots)

exe = qc.compile(prog)<jupyter_output><empty_output><jupyter_text>And now let's estimate $X$ via sampling.<jupyter_code>n_expts = 1000

X_expects = []

for _ in tqdm(range(n_expts)):
    
    samples = qc.run(exe)
    X_eigs = 1 - 2*samples
    X_expect = np.mean(X_eigs)
    X_expects.append(X_expect)
    
plt.figure(figsize=(12, 8))
plt.hist(X_expects, bins=50)
plt.show()<jupyter_output>100%|██████████| 1000/1000 [00:06<00:00, 144.00it/s]<jupyter_text>Similarly, note that $U^{\dagger}ZU=-Y$ where $U = RX(-\pi/2)$, so we can just apply $RX(-\pi/2)$ before measuring, and then carrying out the usual procedure (but taking care to take care of the minus sign).<jupyter_code># ideal value
prog = Program(H(0))
Y_expectation = wfn_sim.expectation(prog, sY(0)+0.0)
print(Y_expectation)
# via sampling
n_shots = 1000

prog = Program()
ro = prog.declare('ro', 'BIT', 1)
prog += H(0)

# rotate into Y-basis
prog += RX(-np.pi/2, 0)

# now measure as usual
prog += MEASURE(0, ro[0])
prog.wrap_in_numshots_loop(n_shots)

exe = qc.compile(prog)

print(exe.program)
n_expts = 1000

Y_expects = []

for _ in tqdm(range(n_expts)):
    
    samples = qc.run(exe)
    Y_eigs = -1.0 * (1 - 2*samples)
    Y_expect = np.mean(Y_eigs)
    Y_expects.append(Y_expect)
    
plt.figure(figsize=(12, 8))
plt.hist(Y_expects, bins=50)
plt.show()<jupyter_output>100%|██████████| 1000/1000 [00:17<00:00, 57.63it/s]<jupyter_text>Readout Error There is a simple way to introduce readout error in pyQuil. The syntax is given [here](https://pyquil-docs.rigetti.com/en/stable/apidocs/autogen/pyquil.quil.Program.define_noisy_readout.html).<jupyter_code>n_shots = 1000

p00 = 0.85
p11 = 0.92

prog = Program()
ro = prog.declare('ro', 'BIT', 1)
prog.define_noisy_readout(0, p00, p11)
prog += H(0)
prog += MEASURE(0, ro[0])
prog.wrap_in_numshots_loop(n_shots)

exe = qc.compile(prog)
n_expts = 1000

Z_expects = []

for _ in tqdm(range(n_expts)):
    
    samples = qc.run(exe)
    Z_eigs = 1 - 2*samples
    Z_expect = np.mean(Z_eigs)
    Z_expects.append(Z_expect)
    
plt.figure(figsize=(12, 8))
plt.hist(Z_expects, bins=50)
plt.show()<jupyter_output>100%|██████████| 1000/1000 [00:16<00:00, 59.24it/s]<jupyter_text>Notice how readout error has caused the center of the Bull curve to shift!  Confusion matrix inversion Let's now repeat the above analysis to measure the observable $ZZ$ in the Bell state.<jupyter_code>n_shots = 1000

# readout noise on qubit 0
p00_q0 = 0.87
p11_q0 = 0.91

# readout noise on qubit 1
p00_q1 = 0.93
p11_q1 = 0.86

noisy_prog = Program()
noisy_prog.define_noisy_readout(0, p00_q0, p11_q0)
noisy_prog.define_noisy_readout(1, p00_q1, p11_q1)

prog = Program()
ro = prog.declare('ro', 'BIT', 2)
prog += noisy_prog
prog += H(0)
prog += CNOT(0,1)
for idx, q in enumerate([0,1]):
    prog += MEASURE(q, ro[idx])
prog.wrap_in_numshots_loop(n_shots)

exe = qc.compile(prog)<jupyter_output><empty_output><jupyter_text>Ideally, since $ZZ \left( \frac{\vert 00 \rangle + \vert 11 \rangle}{\sqrt{2}}\right) = 1$, in the absence of errors the expectation value would be 1. However, readout errors can shift the measured expectation value away from the ideal value, as we observe below.<jupyter_code>qc = get_qc(f'{n_qubits}q-qvm')
n_expts = 1000

ZZ_expects = []

for _ in tqdm(range(n_expts)):
    
    samples = qc.run(exe)
    ZZ_eigs = np.prod(1 - 2*samples, axis=1)
    ZZ_expect = np.mean(ZZ_eigs)
    ZZ_expects.append(ZZ_expect)
    
plt.figure(figsize=(12, 8))
plt.hist(ZZ_expects, bins=50)
plt.show()<jupyter_output>100%|██████████| 1000/1000 [00:33<00:00, 30.28it/s]<jupyter_text>Let us now try to correct for this readout error. As explained in the slides, we can make use of Bayes identity $p(read=x) = \sum_{y} p(read=x \vert sampled=y) \cdot p(sampled=y)$ to infer the true sampling probabilities $p(sampled=y)$. First, we estimate what the empirical probabilities are, i.e. $p(read=x)$.<jupyter_code>qc = get_qc(f'{n_qubits}q-qvm')
samples = qc.run(exe)
samples_str = [''.join(str(q) for q in shot) for shot in samples]

all_bitstrings = list(''.join(bb) for bb in itertools.product(['0','1'], repeat=2))

probs_empirical = []
for bitstring in all_bitstrings:
    frac = samples_str.count(bitstring) / n_shots
    probs_empirical.append(frac)<jupyter_output><empty_output><jupyter_text>Next, we estimate the confusion matrix $p(read=x \vert sampled=y)$. We do this by successively preparing all possible bitstrings $y$, then read them out and measure the fraction of times we obtain each possible bitstring.<jupyter_code>qc = get_qc(f'{n_qubits}q-qvm')
n_expts = 1000

dim = 2**n_qubits
confusion_matrix = np.zeros((dim, dim))

# for (bb1,bb2) in itertools.product(all_bitstrings, repeat=2):
for idx, bb in enumerate(all_bitstrings):
    
    prog = Program()
    prog += noisy_prog
    
    for q, x in enumerate(bb):
        if x=='0':
            pass
        elif x=='1':
            prog += X(q)
            
    for i, q in enumerate([0,1]):
        prog += MEASURE(q, ro[i])

    prog.wrap_in_numshots_loop(n_shots)
    
    exe = qc.compile(prog)
    samples = qc.run(exe)
    samples_str = [''.join(str(q) for q in shot) for shot in samples]
    measured_fracs = [samples_str.count(bitstring)/n_shots for bitstring in all_bitstrings]
    confusion_matrix[:,idx] = measured_fracs<jupyter_output>/Users/sohaib/Documents/pyquil/pyquil/quil.py:1180: UserWarning: Please DECLARE all memory. I'm adding a declaration for the `ro` register, but I won't do this for you in the future.
  "Please DECLARE all memory. I'm adding a declaration for the `ro` register, "<jupyter_text>In vector notation, we can express $p(read=x) = \sum_{y} p(read=x \vert sampled=y) \cdot p(sampled=y)$ as $\vec{p} = P \vec{q}$, where $\vec{p}$ are the empirical probabilities, $P$ is the confusion matrix, and $q$ are the true sampling probabilities. To infer $\vec{q}$, we can simply invert this equation and solve $\vec{q} = P^{-1} \vec{p}$.<jupyter_code>conf_mat_inv = np.linalg.inv(confusion_matrix)
probs_correct = conf_mat_inv @ probs_empirical<jupyter_output><empty_output><jupyter_text>With the corrected probabilities ($\vec{q}$) we can infer the expectation value of the observable $ZZ$.<jupyter_code>ZZ_empirical = (probs_empirical[0] + probs_empirical[3]) - (probs_empirical[1] + probs_empirical[2])
ZZ_corrected = (probs_correct[0] + probs_correct[3]) - (probs_correct[1] + probs_correct[2])

print(f"Empirical ZZ estimate: {ZZ_empirical}")
print(f"Corrected ZZ estimate: {ZZ_corrected}")<jupyter_output>Empirical ZZ estimate: 0.562
Corrected ZZ estimate: 0.8965972964990718<jupyter_text>We see that the corrected estimate is much closer to the ideal value of 1.0. Note however that the corrected "probabilities" are really only quasi-probabilities, as even though they sum to 1<jupyter_code>print(np.sum(probs_correct))<jupyter_output>1.0<jupyter_text>they need not necessarily all be non-negative<jupyter_code>print(probs_correct)<jupyter_output>[0.48376423 0.02600961 0.02569175 0.46453442]
{'fraction_non_alphanumeric': 0.13290208762128786, 'fraction_numerical': 0.04243849848083897, 'mean_word_length': 4.936009307737057, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 63, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '571823', 'license': 'MIT', 'n_tokens_mistral': 4277, 'n_tokens_neox': 3818, 'n_words': 1015, 'path': 'Qunova Computing/rigetti_resources/Sampling and Readout Errors in pyQuil.ipynb', 'repo_name': 'ushnishray/Hackathon2021', 'size': '10203'}