Added documentation to tomography.py

This commit is contained in:
cjwood 2017-08-02 12:33:05 -04:00
parent e503758feb
commit 89cce4d0e7
1 changed files with 368 additions and 149 deletions

View File

@ -2,25 +2,44 @@
# Copyright 2017 IBM RESEARCH. All Rights Reserved.
#
# This file is intended only for use during the USEQIP Summer School 2017.
# Do not distribute.
# It is provided without warranty or conditions of any kind, either express or
# implied.
# An open source version of this file will be included in QISKIT-DEV-PY
# reposity in the future. Keep an eye on the Github repository for updates!
# https://github.com/IBM/qiskit-sdk-py
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# =============================================================================
"""
Functions for performing quantum state and process tomography experiments.
Quantum Tomography Module
Author: Christopher J. Wood <cjwood@us.ibm.com>
Jay Gambetta
TODO:
- Find and fix bug with 2-qubit process tomography basis
- SDP fitting
Description:
This module contains functions for performing quantum state and quantum
process tomography. This includes:
- Functions for generating a set of circuits in a QuantumProgram to
extract tomographically complete sets of measurement data.
- Functions for generating a tomography data set from the QuantumProgram
results after the circuits have been executed on a backend.
- Functions for reconstructing a quantum state, or quantum process
(Choi-matrix) from tomography data sets.
Reconstruction Methods:
Currently implemented reconstruction methods are
- Linear inversion by weighted least-squares fitting.
- Fast maximum likelihood reconstruction using ref [1].
References:
[1] J Smolin, JM Gambetta, G Smith, Phys. Rev. Lett. 108, 070502 (2012).
Open access: arXiv:1106.5458 [quant-ph].
"""
import numpy as np
from functools import reduce
from re import match
@ -41,18 +60,73 @@ where X0 is the projector onto the 0 outcome state of 'X'
"""
def build_state_tomography_circuits(Q_program, name, qubits, qreg, creg):
def build_state_tomography_circuits(Q_program, name, qubits, qreg, creg,
silent=False):
"""
Add state tomography measurement circuits to a QuantumProgram.
The quantum program must contain a circuit 'name', which is treated as a
state preparation circuit. This function then appends the circuit with a
tomographically overcomplete set of measurements in the Pauli basis for
each qubit to be measured. For n-qubit tomography this result in 3 ** n
measurement circuits being added to the quantum program.
Args:
Q_program (QuantumProgram): A quantum program to store the circuits.
name (string): The name of the base circuit to be appended.
qubits (list[int]): a list of the qubit indexes of qreg to be measured.
qreg (QuantumRegister): the quantum register containing qubits to be
measured.
creg (ClassicalRegister): the classical register containing bits to
store measurement outcomes.
silent (bool, optional): hide verbose output.
Returns:
A list of names of the added quantum state tomography circuits.
Example: ['circ_measX0', 'circ_measY0', 'circ_measZ0']
"""
labels = __add_meas_circuits(Q_program, name, qubits, qreg, creg)
print('>> created state tomography circuits for "%s"' % name)
if not silent:
print('>> created state tomography circuits for "%s"' % name)
return labels
# Make private for now, since fit method not yet implemented
def build_process_tomography_circuits(Q_program, name, qubits, qreg, creg):
def build_process_tomography_circuits(Q_program, name, qubits, qreg, creg,
silent=False):
"""
Add process tomography measurement circuits to a QuantumProgram.
The quantum program must contain a circuit 'name', which is the circuit
that will be reconstructed via tomographic measurements. This function
then prepends and appends the circuit with a tomographically overcomplete
set of preparations and measurements in the Pauli basis for
each qubit to be measured. For n-qubit process tomography this result in
(6 ** n) * (3 ** n) circuits being added to the quantum program:
- 3 ** n measurements in the Pauli X, Y, Z bases.
- 6 ** n preparations in the +1 and -1 eigenstates of X, Y, Z.
Args:
Q_program (QuantumProgram): A quantum program to store the circuits.
name (string): The name of the base circuit to be appended.
qubits (list[int]): a list of the qubit indexes of qreg to be measured.
qreg (QuantumRegister): the quantum register containing qubits to be
measured.
creg (ClassicalRegister): the classical register containing bits to
store measurement outcomes.
silent (bool, optional): hide verbose output.
Returns:
A list of names of the added quantum process tomography circuits.
Example:
['circ_prepXp0_measX0', 'circ_prepXp0_measY0', 'circ_prepXp0_measZ0',
'circ_prepXm0_measX0', 'circ_prepXm0_measY0', 'circ_prepXm0_measZ0',
'circ_prepYp0_measX0', 'circ_prepYp0_measY0', 'circ_prepYp0_measZ0',
'circ_prepYm0_measX0', 'circ_prepYm0_measY0', 'circ_prepYm0_measZ0',
'circ_prepZp0_measX0', 'circ_prepZp0_measY0', 'circ_prepZp0_measZ0',
'circ_prepZm0_measX0', 'circ_prepZm0_measY0', 'circ_prepZm0_measZ0']
"""
# add preparation circuits
preps = __add_prep_circuits(Q_program, name, qubits, qreg, creg)
# add measurement circuits for each prep circuit
@ -61,7 +135,8 @@ def build_process_tomography_circuits(Q_program, name, qubits, qreg, creg):
labels += __add_meas_circuits(Q_program, circ, qubits, qreg, creg)
# delete temp prep output
del Q_program._QuantumProgram__quantum_program[circ]
print('>> created process tomography circuits for "%s"' % name)
if not silent:
print('>> created process tomography circuits for "%s"' % name)
return labels
@ -77,31 +152,16 @@ def __tomo_dicts(qubits, basis=None, states=False):
the default is ['X', 'Y', 'Z']
Returns:
a new list of tomo_dict
A new list of tomo_dict
"""
"""
if basis is None:
basis = ['X', 'Y', 'Z']
elif isinstance(basis, dict):
basis = basis.keys()
if not tomos:
return [{qubit: b} for b in basis]
else:
ret = []
for t in tomos:
for b in basis:
t[qubit] = b
ret.append(t.copy())
return ret
"""
if isinstance(qubits, int):
qubits = [qubits]
if basis is None:
basis = __default_basis
basis = __DEFAULT_BASIS
if states is True:
if states:
ns = len(list(basis.values())[0])
lst = [(b, s) for b in basis.keys() for s in range(ns)]
else:
@ -110,14 +170,14 @@ def __tomo_dicts(qubits, basis=None, states=False):
return [dict(zip(qubits, b)) for b in product(lst, repeat=len(qubits))]
def __add_meas_circuits(Q_program, name, qubits, qreg, creg, prep=None):
"""Helper function.
def __add_meas_circuits(Q_program, name, qubits, qreg, creg):
"""
Add measurement circuits to a quantum program.
See: build_state_tomography_circuits.
build_process_tomography_circuits.
"""
if isinstance(qreg, str):
qreg = Q_program.get_quantum_registers(qreg)
if isinstance(creg, str):
creg = Q_program.get_classical_registers(creg)
orig = Q_program.get_circuit(name)
labels = []
@ -128,22 +188,18 @@ def __add_meas_circuits(Q_program, name, qubits, qreg, creg, prep=None):
label = '_meas'
for qubit, op in dic.items():
label += op + str(qubit)
circ = Q_program.create_circuit(label, [qreg], [creg])
circuit = Q_program.create_circuit(label, [qreg], [creg])
# add gates to circuit
for qubit, op in dic.items():
circ.barrier(qreg[qubit])
circuit.barrier(qreg[qubit])
if op == "X":
circ.h(qreg[qubit])
circuit.u2(0., np.pi, qreg[qubit]) # H
elif op == "Y":
circ.sdg(qreg[qubit])
circ.h(qreg[qubit])
if prep:
circ.barrier(qreg[qubit])
else:
circ.measure(qreg[qubit], creg[qubit])
circuit.u2(0., 0.5 * np.pi, qreg[qubit]) # H.S^*
circuit.measure(qreg[qubit], creg[qubit])
# add circuit to QuantumProgram
Q_program.add_circuit(name+label, orig + circ)
Q_program.add_circuit(name+label, orig + circuit)
# add label to output
labels.append(name+label)
# delete temp circuit
@ -153,30 +209,31 @@ def __add_meas_circuits(Q_program, name, qubits, qreg, creg, prep=None):
def __add_prep_gates(circuit, qreg, qubit, op):
"""helper function
"""
Add state preparation gates to a circuit.
"""
p, s = op
if p == "X":
if s == 1:
circuit.x(qreg[qubit])
circuit.h(qreg[qubit])
circuit.u2(np.pi, np.pi, qreg[qubit]) # H.X
else:
circuit.u2(0., np.pi, qreg[qubit]) # H
if p == "Y":
if s == 1:
circuit.x(qreg[qubit])
circuit.h(qreg[qubit])
circuit.s(qreg[qubit])
circuit.u2(-0.5 * np.pi, np.pi, qreg[qubit]) # S.H.X
else:
circuit.u2(0.5 * np.pi, np.pi, qreg[qubit]) # S.H
if p == "Z" and s == 1:
circuit.x(qreg[qubit])
circuit.u3(np.pi, 0., np.pi, qreg[qubit]) # X
def __add_prep_circuits(Q_program, name, qubits, qreg, creg):
"""Helper function.
"""
Add preparation circuits to a quantum program.
See: build_process_tomography_circuits.
"""
if isinstance(qreg, str):
qreg = Q_program.get_quantum_registers(qreg)
if isinstance(creg, str):
creg = Q_program.get_classical_registers(creg)
orig = Q_program.get_circuit(name)
labels = []
@ -206,8 +263,6 @@ def __add_prep_circuits(Q_program, name, qubits, qreg, creg):
###############################################################
# Tomography circuit labels
###############################################################
def __tomo_labels(name, qubits, basis=None, states=False):
"""Helper function.
"""
@ -215,7 +270,7 @@ def __tomo_labels(name, qubits, basis=None, states=False):
state = {0: 'p', 1: 'm'}
for dic in __tomo_dicts(qubits, states=states):
label = ''
if states is True:
if states:
for qubit, op in dic.items():
label += op[0] + state[op[1]] + str(qubit)
else:
@ -225,79 +280,75 @@ def __tomo_labels(name, qubits, basis=None, states=False):
return labels
def state_tomography_labels(name, qubits, basis=None):
def state_tomography_circuit_names(name, qubits):
"""
Return a list of state tomography circuit names.
This list is the same as that returned by the
build_state_tomography_circuits function.
Args:
name (string): the name of the original state preparation
circuit.
qubits: (list[int]): the qubits being measured.
Returns:
A list of circuit names.
"""
return __tomo_labels(name + '_meas', qubits, basis)
return __tomo_labels(name + '_meas', qubits)
# Make private for now, since fit method not yet implemented
def process_tomography_labels(name, qubits, basis=None):
def process_tomography_circuit_names(name, qubits):
"""
Return a list of process tomography circuit names.
This list is the same as that returned by the
build_process_tomography_circuits function.
Args:
name (string): the name of the original circuit to be
reconstructed.
qubits: (list[int]): the qubits being measured.
Returns:
A list of circuit names.
"""
preps = __tomo_labels(name + '_prep', qubits, basis, states=True)
preps = __tomo_labels(name + '_prep', qubits, states=True)
return reduce(lambda acc, c:
acc + __tomo_labels(c + '_meas', qubits, basis),
acc + __tomo_labels(c + '_meas', qubits),
preps, [])
###############################################################
# Tomography measurement outcomes
# Preformatting count data
###############################################################
# Note: So far only for state tomography
# Default Pauli basis
__default_basis = {'X': [np.array([[0.5, 0.5],
[0.5, 0.5]]),
np.array([[0.5, -0.5],
[-0.5, 0.5]])],
'Y': [np.array([[0.5, -0.5j],
[0.5j, 0.5]]),
np.array([[0.5, 0.5j],
[-0.5j, 0.5]])],
'Z': [np.array([[1, 0],
[0, 0]]),
np.array([[0, 0],
[0, 1]])]}
def __counts_keys(n):
"""helper function.
"""Generate outcome bitstrings for n-qubits.
Args:
n (int): the number of qubits.
Returns:
A list of bitstrings ordered as follows:
Example: n=2 returns ['00', '01', '10', '11'].
"""
return [bin(j)[2:].zfill(n) for j in range(2 ** n)]
def __get_meas_basis_ops(tup, basis):
# reverse tuple so least significant qubit is to the right
return reduce(lambda acc, b: [np.kron(a, j)
for a in acc for j in basis[b]],
reversed(tup), [1])
def __meas_basis(n, basis):
return [dict(zip(__counts_keys(n), __get_meas_basis_ops(key, basis)))
for key in product(basis.keys(), repeat=n)]
def __get_prep_basis_op(dic, basis):
keys = sorted(dic.keys()) # order qubits [0,1,...]
tups = [dic[k] for k in keys]
return reduce(lambda acc, b: np.kron(basis[b[0]][b[1]], acc),
tups, [1])
def __prep_basis(n, basis):
# use same function as prep circuits to get order right
ordered = __tomo_dicts(range(n), states=True)
return [__get_prep_basis_op(dic, basis) for dic in ordered]
def marginal_counts(counts, meas_qubits):
"""
Returns a list of the measurement outcome strings for meas_qubits qubits
in an n-qubit system.
Compute the marginal counts for a subset of measured qubits.
Args:
counts (dict{str:int}): the counts returned from a backend.
meas_qubits (list[int]): the qubits to return the marginal
counts distribution for.
Returns:
A counts dict for the meas_qubits.abs
Example: if counts = {'00': 10, '01': 5}
marginal_counts(counts, [0]) returns {'0': 15, '1': 0}.
marginal_counts(counts, [0]) returns {'0': 10, '1': 5}.
"""
# Extract total number of qubits from count keys
@ -326,13 +377,97 @@ def marginal_counts(counts, meas_qubits):
return dict(zip(meas_keys, meas_counts))
def state_tomography_data(Q_program, name, meas_qubits,
backend=None, basis=None):
###############################################################
# Tomography preparation and measurement bases
###############################################################
# Default Pauli basis
# This corresponds to measurements in the X, Y, Z basis where
# Outcomes 0,1 are the +1,-1 eigenstates respectively.
# State preparation is also done in the +1 and -1 eigenstates.
__DEFAULT_BASIS = {'X': [np.array([[0.5, 0.5],
[0.5, 0.5]]),
np.array([[0.5, -0.5],
[-0.5, 0.5]])],
'Y': [np.array([[0.5, -0.5j],
[0.5j, 0.5]]),
np.array([[0.5, 0.5j],
[-0.5j, 0.5]])],
'Z': [np.array([[1, 0],
[0, 0]]),
np.array([[0, 0],
[0, 1]])]}
def __get_meas_basis_ops(tup, basis):
"""
Return a n-qubit projector for a given measurement.
"""
# reverse tuple so least significant qubit is to the right
return reduce(lambda acc, b: [np.kron(a, j)
for a in acc for j in basis[b]],
reversed(tup), [1])
def __meas_basis(n, basis):
"""
Return an ordered list of n-qubit measurment projectors.
"""
return [dict(zip(__counts_keys(n), __get_meas_basis_ops(key, basis)))
for key in product(basis.keys(), repeat=n)]
def __get_prep_basis_op(dic, basis):
"""
Return an n-qubit projector for a given prepration.
"""
keys = sorted(dic.keys()) # order qubits [0,1,...]
tups = [dic[k] for k in keys]
return reduce(lambda acc, b: np.kron(basis[b[0]][b[1]], acc),
tups, [1])
def __prep_basis(n, basis):
"""
Return an ordered list of n-qubit preparation projectors.
"""
# use same function as prep circuits to get order right
ordered = __tomo_dicts(range(n), states=True)
return [__get_prep_basis_op(dic, basis) for dic in ordered]
def state_tomography_data(Q_program, name, meas_qubits, backend=None,
basis=None):
"""
Return a list of state tomography measurement outcomes.
Args:
Q_program (QuantumProgram): A quantum program containing count data
from execution on a backend.
name (string): The name of the base state preparation circuit.
meas_qubits (list[int]): a list of the qubit indexes measured.
backend (str, optional): the backend the tomography circuits were
executed on. If not specified it will use the
default backend of QuantumProgram.get_counts.
basis (basis dict, optional): the basis used for measurement. Default
is the Pauli basis.
Returns:
A list of dicts for the outcome of each state tomography
measurement circuit. The keys of the dictionary are
{
'counts': dict('str': int),
<the marginal counts for measured qubits>,
'shots': int,
<total number of shots for measurement circuit>
'meas_basis': dict('str': np.array)
<the projector for the measurement outcomes>
}
"""
if basis is None:
basis = __default_basis
labels = state_tomography_labels(name, meas_qubits)
basis = __DEFAULT_BASIS
labels = state_tomography_circuit_names(name, meas_qubits)
counts = [marginal_counts(Q_program.get_counts(circ, backend), meas_qubits)
for circ in labels]
shots = [sum(c.values()) for c in counts]
@ -342,14 +477,40 @@ def state_tomography_data(Q_program, name, meas_qubits,
return ret
def process_tomography_data(Q_program, name, meas_qubits,
backend=None, basis=None):
def process_tomography_data(Q_program, name, meas_qubits, backend=None,
basis=None):
"""
Return a list of process tomography measurement outcomes.
Args:
Q_program (QuantumProgram): A quantum program containing count data
from execution on a backend.
name (string): The name of the circuit being reconstructed.
meas_qubits (list[int]): a list of the qubit indexes measured.
backend (str, optional): the backend the tomography circuits were
executed on. If not specified it will use the
default backend of QuantumProgram.get_counts.
basis (basis dict, optional): the basis used for measurement. Default
is the Pauli basis.
Returns:
A list of dicts for the outcome of each process tomography
measurement circuit. The keys of the dictionary are
{
'counts': dict('str': int),
<the marginal counts for measured qubits>,
'shots': int,
<total number of shots for measurement circuit>
'meas_basis': dict('str': np.array),
<the projector for the measurement outcomes>
'prep_basis': np.array,
<the projector for the prepared input state>
}
"""
if basis is None:
basis = __default_basis
basis = __DEFAULT_BASIS
n = len(meas_qubits)
labels = process_tomography_labels(name, meas_qubits)
labels = process_tomography_circuit_names(name, meas_qubits)
counts = [marginal_counts(Q_program.get_counts(circ, backend), meas_qubits)
for circ in labels]
shots = [sum(c.values()) for c in counts]
@ -364,15 +525,17 @@ def process_tomography_data(Q_program, name, meas_qubits,
dic['shots'] = sts
return ret
###############################################################
# Tomographic Reconstruction functions.
###############################################################
def __tomo_basis_matrix(meas_basis):
"""Return a matrix of vectorized measurement operators.
Measurement operators S = sum_j |j><M_j|.
Args:
meas_basis(list(array_like)): measurement operators [M_j].
Returns:
The operators S = sum_j |j><M_j|.
"""
n = len(meas_basis)
d = meas_basis[0].size
@ -382,6 +545,17 @@ def __tomo_basis_matrix(meas_basis):
def __tomo_linear_inv(freqs, ops, weights=None, trace=None):
"""
Reconstruct a matrix through linear inversion.
Args:
freqs (list[float]): list of observed frequences.
ops (list[np.array]): list of corresponding projectors.
weights (list[float] or array_like, optional):
weights to be used for weighted fitting.
trace (float, optional): trace of returned operator.
Returns:
A numpy array of the reconstructed operator.
"""
# get weights matrix
if weights is not None:
@ -410,16 +584,33 @@ def __tomo_linear_inv(freqs, ops, weights=None, trace=None):
return ret
def __leastsq_fit(data, weights=None, trace=None, beta=0.5):
def __leastsq_fit(data, weights=None, trace=None, beta=None):
"""
Reconstruct a state from unconstrained least-squares fitting.
Args:
data (list[dict]): state or process tomography data.
weights (list or array, optional): weights to use for least squares
fitting. The default is standard deviation from a binomial
distribution.
trace (float, optional): trace of returned operator. The default is 1.
beta (float >=0, optional): hedge parameter for computing frequencies
from zero-count data. The default value is 0.50922.
Returns:
A numpy array of the reconstructed operator.
"""
if trace is None:
trace = 1. # default to unit trace
ks = data[0]['counts'].keys()
K = len(ks)
# Get counts and shots
ns = np.array([dat['counts'][k] for dat in data for k in ks])
shots = np.array([dat['shots'] for dat in data for k in ks])
# convert to freqs using beta to hedge against zero counts
if beta is None:
beta = 0.50922
freqs = (ns + beta) / (shots + K * beta)
# Use standard least squares fitting weights
@ -438,13 +629,25 @@ def __leastsq_fit(data, weights=None, trace=None, beta=0.5):
return __tomo_linear_inv(freqs, ops, weights, trace=trace)
def __wizard(rho, epsilon=0.):
def __wizard(rho, epsilon=None):
"""
Maps an operator to the nearest positive semidefinite operator
Returns the nearest postitive semidefinite operator to an operator.
This method is based on reference [1]. It constrains positivity
by setting negative eigenvalues to zero and rescaling the positive
eigenvalues.
See arXiv:1106.5458 [quant-ph]
Args:
rho (array_like): the input operator.
epsilon(float >=0, optional): threshold for truncating small
eigenvalues values to zero.
Returns:
A positive semidefinite numpy array.
"""
if epsilon is None:
epsilon = 0. # default value
dim = len(rho)
rho_wizard = np.zeros([dim, dim])
v, w = np.linalg.eigh(rho) # v eigenvecrors v[0] < v[1] <...
@ -463,15 +666,23 @@ def __wizard(rho, epsilon=0.):
def __get_option(opt, options):
"""
Return an optional value or None if not found.
"""
if options is not None:
if opt in options:
return options[opt]
return None
def fit_tomography_data(data, method='wizard', options=None):
def fit_tomography_data(data, method=None, options=None):
"""
Reconstruct a state or Choi-matrix from tomography data.
Reconstruct a density matrix or process-matrix from tomography data.
If the input data is state_tomography_data the returned operator will
be a density matrix. If the input data is process_tomography_data the
returned operator will be a Choi-matrix in the column-vectorization
convention.
Args:
data (dict): process tomography measurement data.
@ -483,27 +694,35 @@ def fit_tomography_data(data, method='wizard', options=None):
Returns:
The fitted operator.
Available methods:
- 'wizard' (Default): The returned operator will be constrained to be
positive-semidefinite.
Options:
- 'trace': the trace of the returned operator.
The default value is 1.
- 'beta': hedging parameter for computing frequencies from
zero-count data. The default value is 0.50922.
- 'epsilon: threshold for truncating small eigenvalues to zero.
The default value is 0
- 'leastsq': Fitting without postive-semidefinite constraint.
Options:
- 'trace': Same as for 'wizard' method.
- 'beta': Same as for 'wizard' method.
"""
if method is None:
method = 'wizard' # set default method
if method in ['wizard', 'leastsq']:
# get options
trace = __get_option('trace', options)
if trace is None:
trace = 1
beta = __get_option('beta', options)
if beta is None:
beta = 0.50922
# fit state
rho = __leastsq_fit(data, trace=trace, beta=beta)
if method == 'wizard':
# Use wizard method to constrain positivity
epsilon = __get_option('epsilon', options)
if epsilon is None:
epsilon = 0.
rho = __wizard(rho, epsilon=epsilon)
return rho
else:
# TODO: raise exception for unknown method
print('error: method unrecongnized')
pass
return rho
print('error: method unknown reconstruction method "%s"' % method)