This notebook walks through how to run Bell state tomography on a noisy QVM, using parametric compilation and pyQuil's Experiment
framework. This notebook is copied partially from the rigetti/qcs-paper repository, where it was used to produce Figure A2 from A quantum-classical cloud platform optimized for variational hybrid algorithms.
import itertools
from typing import Generator, List
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.cm import ScalarMappable
from matplotlib.colors import LinearSegmentedColormap, Normalize
from scipy.linalg import pinv
import pyquil.simulation.matrices as psm
from forest.benchmarking.distance_measures import fidelity
from forest.benchmarking.operator_tools.superoperator_transformations import unvec, vec
from forest.benchmarking.utils import all_traceless_pauli_terms
from pyquil import get_qc, Program
from pyquil.experiment import (
Experiment, ExperimentResult, ExperimentSetting, correct_experiment_result, zeros_state
)
from pyquil.gates import CNOT, H, RESET
from pyquil.paulis import PauliTerm
from pyquil.simulation.tools import lifted_pauli
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) /tmp/ipykernel_63/3981836547.py in <module> 10 11 import pyquil.simulation.matrices as psm ---> 12 from forest.benchmarking.distance_measures import fidelity 13 from forest.benchmarking.operator_tools.superoperator_transformations import unvec, vec 14 from forest.benchmarking.utils import all_traceless_pauli_terms /usr/local/lib/python3.7/site-packages/forest/benchmarking/distance_measures.py in <module> 5 from scipy.linalg import fractional_matrix_power 6 from scipy.optimize import minimize_scalar ----> 7 from forest.benchmarking.operator_tools.calculational import sqrtm_psd 8 9 /usr/local/lib/python3.7/site-packages/forest/benchmarking/operator_tools/__init__.py in <module> 2 from .channel_approximation import * 3 from .compose_superoperators import * ----> 4 from .project_superoperators import * 5 from .superoperator_transformations import * 6 from .random_operators import * /usr/local/lib/python3.7/site-packages/forest/benchmarking/operator_tools/project_superoperators.py in <module> 14 from scipy import linalg 15 from forest.benchmarking.operator_tools.calculational import partial_trace ---> 16 from forest.benchmarking.operator_tools.superoperator_transformations import vec, unvec, kraus2choi 17 18 /usr/local/lib/python3.7/site-packages/forest/benchmarking/operator_tools/superoperator_transformations.py in <module> 28 from typing import Sequence, Tuple, List 29 import numpy as np ---> 30 from forest.benchmarking.utils import n_qubit_pauli_basis 31 32 /usr/local/lib/python3.7/site-packages/forest/benchmarking/utils.py in <module> 10 11 from pyquil.gates import I, RX, RY, RZ, H, MEASURE ---> 12 from pyquil.gate_matrices import X, Y, Z 13 from pyquil.paulis import PauliTerm 14 from pyquil.quil import Program ModuleNotFoundError: No module named 'pyquil.gate_matrices'
qubits = (0, 1)
shots = 500
qc = get_qc("2q-noisy-qvm")
Experiment
¶p = Program()
p += RESET()
p += H(qubits[0])
p += CNOT(qubits[0], qubits[1])
p.wrap_in_numshots_loop(shots)
print(p)
def state_tomo_settings(qubits: List[int]) -> Generator[ExperimentSetting, None, None]:
"""
Adapted from forest.benchmarking.tomography._state_tomo_settings,
to use pyquil.experiment.ExperimentSetting objects instead.
"""
list_of_terms = all_traceless_pauli_terms(qubits)
for obs in all_traceless_pauli_terms(qubits):
yield ExperimentSetting(
in_state=zeros_state(qubits),
out_operator=obs,
)
state_tomography_settings = list(state_tomo_settings(qubits))
state_tomography_settings
bell_state_tomography = Experiment(settings=state_tomography_settings, program=p)
bell_state_tomography
%%time
results = qc.experiment(bell_state_tomography)
results
%%time
calibrations = qc.calibrate(bell_state_tomography)
calibrations
results_corrected = []
for r, c in zip(results, calibrations):
results_corrected.append(correct_experiment_result(r, c))
results_corrected
def build_rho_true() -> np.ndarray:
"""
Generate the density matrix for state |00> + |11>.
"""
psi00 = np.array([[1], [0], [0], [0]])
bell00 = psm.CNOT @ np.kron(psm.H, psm.I) @ psi00
return np.outer(bell00, bell00.T.conj())
rho_true = build_rho_true()
rho_true
def linear_inv_state_estimate(results: List[ExperimentResult], qubits: List[int]) -> np.ndarray:
"""
Adapted from forest.benchmarking.tomography.linear_inv_state_estimate,
to use pyquil.experiment.ExperimentResult objects instead.
"""
measurement_matrix = np.vstack([
vec(lifted_pauli(result.setting.out_operator, qubits=qubits)).T.conj() for result in results])
expectations = np.array([result.expectation for result in results])
rho = pinv(measurement_matrix) @ expectations
dim = 2**len(qubits)
return unvec(rho) + np.eye(dim) / dim
rho_simulated = linear_inv_state_estimate(results_corrected, [0,1])
np.round(rho_simulated, 3)
DARK_TEAL = "#47717d"
GOLD = "#f8ba2b"
LIGHT_TEAL = "#66acb4"
NAVY = "#00507b"
# build a custom colormap for the hinton plot
lsc = LinearSegmentedColormap.from_list(name="rigetti", colors=[NAVY, GOLD, LIGHT_TEAL, DARK_TEAL])
ANGLE_MAPPER = ScalarMappable(norm=Normalize(vmin=-np.pi, vmax=np.pi))
ANGLE_MAPPER.set_cmap(lsc)
def hinton(matrix: np.ndarray, ax: plt.Axes) -> None:
"""
Adapted from forest.benchmarking.tomography.hinton to use custom colors.
"""
max_weight=1.0
ax.patch.set_facecolor("white")
ax.set_aspect("equal", "box")
ax.xaxis.set_major_locator(plt.NullLocator())
ax.yaxis.set_major_locator(plt.NullLocator())
for (x, y), w in np.ndenumerate(matrix):
color = np.arctan2(w.real, w.imag)
color = ANGLE_MAPPER.to_rgba(color)
size = np.sqrt(np.abs(w) / max_weight)
rect = plt.Rectangle([x - size / 2, y - size / 2], size, size,
facecolor=color, edgecolor=color)
ax.add_patch(rect)
ax.set_xlim((-max_weight / 2, matrix.shape[0] - max_weight / 2))
ax.set_ylim((-max_weight / 2, matrix.shape[1] - max_weight / 2))
ax.autoscale_view()
ax.invert_yaxis()
fig_qvm, (ax0_qvm, ax1_qvm) = plt.subplots(2, 1, figsize=(6, 12))
fig_qvm.subplots_adjust(hspace=0.1)
hinton(rho_true, ax=ax0_qvm)
hinton(rho_simulated, ax=ax1_qvm)
print(f"Simulated Bell state fidelity = {np.round(fidelity(rho_true, rho_simulated), 4)*100:.2f}%")