<jupyter_start><jupyter_text>State tomographyState tomography involves measuring a quantum state repeatedly in the bases given by `itertools.product(['X', 'Y', 'Z], repeat=n_qubits)`. From these measurements, we can reconstruct a density matrix $\rho$.<jupyter_code>import numpy as np
from pyquil import Program, get_qc
from pyquil.gates import *<jupyter_output><empty_output><jupyter_text>Construct a state with a `Program`We'll construct a two-qubit graph state by Hadamarding all qubits and then applying a controlled-Z operation across edges of our graph. In the two-qubit case, there's only one edge.<jupyter_code>qubits = [0, 1]
program = Program()
for qubit in qubits:
    program += H(qubit)
program += CZ(qubits[0], qubits[1])
print(program)<jupyter_output><empty_output><jupyter_text>Construct a `TomographyExperiment` for state tomographyWe can print this out to see the 16 measurements we will perform.<jupyter_code>from forest_benchmarking.tomography import generate_state_tomography_experiment
experiment = generate_state_tomography_experiment(program=program, qubits=qubits)
print(experiment)<jupyter_output><empty_output><jupyter_text>Optional groupingWe can simultaneously estimate some of these observables<jupyter_code>from pyquil.operator_estimation import group_experiments
print(group_experiments(experiment))<jupyter_output><empty_output><jupyter_text>PyQuil will run the tomography programs<jupyter_code>from pyquil.operator_estimation import measure_observables

qc = get_qc('2q-pyqvm')
# Over-write full quilc compilation with a much more simple
# version that *only* substitutes gates to Rigetti-native gates.
# We don't want to accidentally compile away our tomography circuit
# or map to different qubits.
from forest_benchmarking.compilation import basic_compile
qc.compiler.quil_to_native_quil = basic_compile

results = list(measure_observables(qc=qc, tomo_experiment=experiment, n_shots=100_000))
results<jupyter_output><empty_output><jupyter_text>We can look at a bunch of numbers...<jupyter_code>from forest_benchmarking.tomography import linear_inv_state_estimate
rho = linear_inv_state_estimate(results, qubits=qubits)
np.round(rho, 3)
psi = (1/2) * np.array([1, 1, 1, -1])
rho_true = np.outer(psi, psi.T.conj())
rho_true<jupyter_output><empty_output><jupyter_text>Or visualize using Hinton plots<jupyter_code>from matplotlib import pyplot as plt
from forest_benchmarking.plotting import hinton
fig, (ax1, ax2) = plt.subplots(1, 2)
hinton(rho_true, ax=ax1)
hinton(rho, ax=ax2)
ax1.set_title('Analytical')
ax2.set_title('Estimated')
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Matrix norm between true and estimated is low<jupyter_code>np.linalg.norm(rho - rho_true)<jupyter_output><empty_output><jupyter_text>Linear inversion estimate<jupyter_code>from forest_benchmarking.tomography import linear_inv_state_estimate
rho = linear_inv_state_estimate(results, qubits=qubits)

print(np.round(rho, 4))
print('Purity =', np.trace(rho @ rho))
hinton(rho)<jupyter_output><empty_output><jupyter_text>Maximum Liklihood Estimate (MLE) via diluted iterative method<jupyter_code>from forest_benchmarking.tomography import iterative_mle_state_estimate
est_mle, status = iterative_mle_state_estimate(results=results, qubits=qubits,
                                               dilution=0.5)
rho = est_mle.estimate.state_point_est
print(np.around(rho, decimals=4))
print('Purity =', np.trace(rho @ rho))
hinton(rho)<jupyter_output><empty_output><jupyter_text>MLE with Max Entropy constraint<jupyter_code>est_max_ent, stat = iterative_mle_state_estimate(results=results, qubits=qubits,
                                                 dilution=0.5, entropy_penalty=0.005)
rho = est_max_ent.estimate.state_point_est
print(np.around(rho, decimals=4))
print('Purity =', np.trace(rho @ rho))
hinton(rho)<jupyter_output><empty_output><jupyter_text>MLE with Hedging parameter<jupyter_code>est_hedging, stat = iterative_mle_state_estimate(results=results, qubits=qubits,
                                                 dilution=0.5, beta=.61)
rho = est_hedging.estimate.state_point_est
print(np.around(rho, decimals=4))
print('Purity = ', np.trace(rho @ rho))
hinton(rho)<jupyter_output><empty_output><jupyter_text>Project an unphysical state to the closest physical state<jupyter_code>from forest_benchmarking.tomography import project_density_matrix
rho_unphys = np.random.uniform(-1, 1, (2, 2)) \
    * np.exp(1.j * np.random.uniform(-np.pi, np.pi, (2, 2)))
rho_phys = project_density_matrix(rho_unphys)

fig, (ax1, ax2) = plt.subplots(1, 2)
hinton(rho_unphys, ax=ax1)
hinton(rho_phys, ax=ax2)
ax1.set_title('Unphysical')
ax2.set_title('Physical projection')
fig.tight_layout()
# Test the wizard method. Example from fig 1 of maximum likelihood minimum effort 
# https://doi.org/10.1103/PhysRevLett.108.070502

eigs = np.diag(np.array(list(reversed([3.0/5, 1.0/2, 7.0/20, 1.0/10, -11.0/20]))))
phys = project_density_matrix(eigs)
np.allclose(phys, np.diag([0, 0, 1.0/5, 7.0/20, 9.0/20]))<jupyter_output><empty_output><jupyter_text>Lightweight Bootstrap for functionals of the state<jupyter_code>import forest_benchmarking.distance_measures as dm
from forest_benchmarking.tomography import estimate_variance
def my_mle_estimator(_results, _qubits):
    return iterative_mle_state_estimate(_results, _qubits, dilution=0.5, entropy_penalty=0.0, beta=0.0)[0]<jupyter_output><empty_output><jupyter_text>**Purity**<jupyter_code>mle_est = estimate_variance(results, qubits, my_mle_estimator, dm.purity,
                            n_resamples=40, project_to_physical=True)
lin_inv_est = estimate_variance(results, qubits, linear_inv_state_estimate, dm.purity,
                                n_resamples=40, project_to_physical=True)
print(mle_est)
print(lin_inv_est)<jupyter_output><empty_output><jupyter_text>**Fidelity**<jupyter_code>mle_est = estimate_variance(results, qubits, my_mle_estimator, dm.fidelity,
                            target_state=rho_true, n_resamples=40, project_to_physical=True)
lin_inv_est = estimate_variance(results, qubits, linear_inv_state_estimate, dm.fidelity,
                                target_state=rho_true, n_resamples=40, project_to_physical=True)
print(mle_est)
print(lin_inv_est)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10308622335348393, 'fraction_numerical': 0.019408208717785556, 'mean_word_length': 6.226436781609196, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 65, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '110397', 'license': 'Apache-2.0', 'n_tokens_mistral': 2265, 'n_tokens_neox': 2129, 'n_words': 422, 'path': 'examples/tomography_state.ipynb', 'repo_name': 'msohaibalam/forest-benchmarking', 'size': '6286'}