scnc
/
qiskit
mirror of https://github.com/Qiskit/qiskit.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

updated state tomography, added visualizations library, added qi library 

- state tomography supports least squares fitting, and wizard method
- qi contains partial trace, state_fidelity, and in future will contain
other common quantum info functions
This commit is contained in:
cjwood 2017-07-28 12:18:59 -04:00
parent 8a3279791f
commit dd146c7924
5 changed files with 1699 additions and 859 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

241
tools/qi.py Normal file
View File

@ -0,0 +1,241 @@
# -*- coding: utf-8 -*-
# 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
# =============================================================================
"""
A collection of useful quantum information functions.
Author: Christopher J. Wood <cjwood@us.ibm.com>
Currently this file is very sparse. More functions will be added in
the future.
"""
import numpy as np
from scipy.linalg import sqrtm
from tools.pauli import pauli_group
###############################################################
# State manipulation.
###############################################################
def partial_trace(state, sys, dims=None):
"""
Partial trace over subsystems of multi-partite matrix.
Note that subsystems are ordered as rho012 = rho0(x)rho1(x)rho2.
Args:
state (NxN matrix_like): a matrix
sys (list(int): a list of subsystems (starting from 0) to trace over
dims (list(int), optional): a list of the dimensions of the subsystems.
If this is not set it will assume all subsystems are qubits.
Returns:
A matrix with the appropriate subsytems traced over.
"""
# convert op to density matrix
rho = np.array(state)
if rho.ndim == 1:
rho = outer(rho) # convert state vector to density mat
# compute dims if not specified
if dims is None:
n = int(np.log2(len(rho)))
dims = [2 for i in range(n)]
if len(rho) != 2 ** n:
print("ERRROR!")
else:
dims = list(dims)
# reverse sort trace sys
if isinstance(sys, int):
sys = [sys]
else:
sys = sorted(sys, reverse=True)
# trace out subsystems
for j in sys:
# get trace dims
dpre = dims[:j]
dpost = dims[j+1:]
dim1 = int(np.prod(dpre))
dim2 = int(dims[j])
dim3 = int(np.prod(dpost))
# dims with sys-j removed
dims = dpre + dpost
# do the trace over j
rho = __trace_middle(rho, dim1, dim2, dim3)
return rho
def __trace_middle(op, dim1=1, dim2=1, dim3=1):
"""
Partial trace over middle system of tripartite state.
Args:
op (NxN matrix_like): a tri-partite matrix
dim1: dimension of the first subsystem
dim2: dimension of the second (traced over) subsystem
dim3: dimension of the third subsystem
Returns:
A (D,D) matrix where D = dim1 * dim3
"""
op = op.reshape(dim1, dim2, dim3, dim1, dim2, dim3)
d = dim1 * dim3
return op.trace(axis1=1, axis2=4).reshape(d, d)
def vectorize(rho, basis='col'):
"""Flatten an operator to a vector in a specified basis.
Args:
rho (ndarray): a density matrix.
basis (str): the basis to vectorize in. Allowed values are
- 'col' (default) flattens to column-major vector.
- 'row' flattens to row-major vector.
- 'pauli'flattens in the n-qubit Pauli basis.
Returns:
ndarray: the resulting vector.
"""
if basis == 'col':
return rho.flatten(order='F')
elif basis == 'row':
return rho.flatten(order='C')
elif basis == 'pauli':
num = int(np.log2(len(rho))) # number of qubits
if len(rho) != 2**num:
print('error input state must be n-qubit state')
vals = list(map(lambda x: np.real(np.trace(np.dot(x.to_matrix(), rho))),
pauli_group(num)))
return np.array(vals)
def devectorize(vec, basis='col'):
"""Devectorize a vectorized square matrix.
Args:
vec (ndarray): a vectorized density matrix.
basis (str): the basis to vectorize in. Allowed values are
- 'col' (default): flattens to column-major vector.
- 'row': flattens to row-major vector.
- 'pauli': flattens in the n-qubit Pauli basis.
Returns:
ndarray: the resulting matrix.
"""
d = int(np.sqrt(vec.size)) # the dimension of the matrix
if len(vec) != d*d:
print('error input is not a vectorized square matrix')
if basis == 'col':
return vec.reshape(d, d, order='F')
elif basis == 'row':
return vec.reshape(d, d, order='C')
elif basis == 'pauli':
num = int(np.log2(d)) # number of qubits
if d != 2 ** num:
print('error input state must be n-qubit state')
pbasis = np.array([p.to_matrix() for p in pauli_group(2)]) / num**2
return np.tensordot(vec, pbasis, axes=1)
def chop(op, epsilon=1e-10):
"""
Truncate small values of a complex array.
Args:
op (array_like): array to truncte small values.
epsilon (float): threshold.
Returns:
A new operator with small values set to zero.
"""
op.real[abs(op.real) < epsilon] = 0.0
op.imag[abs(op.imag) < epsilon] = 0.0
return op
def outer(v1, v2=None):
"""
Construct the outer product of two vectors.
The second vector argument is optional, if absent the projector
of the first vector will be returned.
Args:
v1 (ndarray): the first vector.
v2 (ndarray): the (optional) second vector.
Returns:
The matrix |v1><v2|.
"""
if v2 is None:
u = np.array(v1).conj()
else:
u = np.array(v2).conj()
return np.outer(v1, u)
###############################################################
# Measures.
###############################################################
def state_fidelity(state1, state2):
"""Return the state fidelity between two quantum states.
Either input may be a state vector, or a density matrix.
Args:
state1: a quantum state vector or density matrix.
state2: a quantum state vector or density matrix.
Returns:
The state fidelity F(state1, state2).
"""
# convert input to numpy arrays
s1 = np.array(state1)
s2 = np.array(state2)
# fidelity of two state vectors
if s1.ndim == 1 and s2.ndim == 1:
return np.abs(np.dot(s2.conj(), s1))
# fidelity of vector and density matrix
elif s1.ndim == 1:
# psi = s1, rho = s2
return np.sqrt(np.abs(np.dot(s1.conj(), np.dot(s2, s1))))
elif s2.ndim == 1:
# psi = s2, rho = s1
return np.sqrt(np.abs(np.dot(s2.conj(), np.dot(s1, s2))))
# fidelity of two density matrices
else:
(s1sq, err1) = sqrtm(s1, disp=False)
(s2sq, err2) = sqrtm(s2, disp=False)
return np.linalg.norm(np.dot(s1sq, s2sq), ord='nuc')
###############################################################
# Other.
###############################################################
def is_pos_def(x):
return np.all(np.linalg.eigvals(x) > 0)

848
tools/tomography.py
View File

@ -17,468 +17,280 @@ from Smolin, Gambetta, Smith Phys. Rev. Lett. 108, 070502 (arXiv: 1106.5458)
Author: Christopher J. Wood <cjwood@us.ibm.com>
Jay Gambetta
Andrew Cross
TODO: Process tomography, SDP fitting
"""
import numpy as np
from functools import reduce
from scipy import linalg as la
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch
from tools.pauli import pauli_group, pauli_singles
from re import match
from itertools import product
from qiskit import QuantumProgram
from tools.qi import vectorize, devectorize, outer
###############################################################
# Plotting tools
# NEW CIRCUIT GENERATION
###############################################################
class Arrow3D(FancyArrowPatch):
"""Standard 3D arrow."""
def __init__(self, xs, ys, zs, *args, **kwargs):
"""Create arrow."""
FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
self._verts3d = xs, ys, zs
def draw(self, renderer):
"""Draw the arrow."""
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
FancyArrowPatch.draw(self, renderer)
"""
Basis should be specified as a dictionary {'X': [X0, X1], 'Y': [Y0, Y1], 'Z': [Z0, Z1]}
where X0 is the projector onto the 0 outcome state of 'X'
"""
# COMPLEMENT = {'1': '0', '0': '1'}
def build_state_tomography_circuits(Q_program, name, qubits, qreg, creg):
"""
"""
labels = __add_meas_circuits(Q_program, name, qubits, qreg, creg)
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):
"""
"""
# add preparation circuits
preps = __add_prep_circuits(Q_program, name, qubits, qreg, creg)
# add measurement circuits for each prep circuit
labels = []
for circ in preps:
labels += __add_meas_circuits(Q_program, circ, qubits, qreg, creg)
# delete temp prep output
del Q_program._QuantumProgram__quantum_program['circuits'][circ]
print('>> created process tomography circuits for "%s"' % name)
return labels
# def compliment(value):
# """Swap 1 and 0 in a vector."""
# return ''.join(COMPLEMENT[x] for x in value)
def __tomo_dicts(qubits, basis=None):
"""Helper function.
def plot_bloch_vector(bloch, title=""):
"""Plot the Bloch sphere.
Plot a sphere, axes, the Bloch vector, and its projections onto each axis.
Build a dictionary assigning a basis element to a qubit.
Args:
bloch (list[double]): array of three elements where [<x>, <y>,<z>]
title (str): a string that represents the plot title
qubit (int): the qubit to add
tomos (list[dict]): list of tomo_dicts to add to
basis (list[str], optional): basis to use. If not specified
the default is ['X', 'Y', 'Z']
Returns:
none: plot is shown with matplotlib to the screen
a new list of tomo_dict
"""
# Set arrow lengths
arlen = 1.3
# Plot semi-transparent sphere
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='3d')
ax.set_aspect("equal")
ax.plot_surface(x, y, z, color=(.5, .5, .5), alpha=0.1)
# Plot arrows (axes, Bloch vector, its projections)
xa = Arrow3D([0, arlen], [0, 0], [0, 0], mutation_scale=20, lw=1,
arrowstyle="-|>", color=(.5, .5, .5))
ya = Arrow3D([0, 0], [0, arlen], [0, 0], mutation_scale=20, lw=1,
arrowstyle="-|>", color=(.5, .5, .5))
za = Arrow3D([0, 0], [0, 0], [0, arlen], mutation_scale=20, lw=1,
arrowstyle="-|>", color=(.5, .5, .5))
a = Arrow3D([0, bloch[0]], [0, bloch[1]], [0, bloch[2]], mutation_scale=20,
lw=2, arrowstyle="simple", color="k")
bax = Arrow3D([0, bloch[0]], [0, 0], [0, 0], mutation_scale=20, lw=2,
arrowstyle="-", color="r")
bay = Arrow3D([0, 0], [0, bloch[1]], [0, 0], mutation_scale=20, lw=2,
arrowstyle="-", color="g")
baz = Arrow3D([0, 0], [0, 0], [0, bloch[2]], mutation_scale=20, lw=2,
arrowstyle="-", color="b")
arrowlist = [xa, ya, za, a, bax, bay, baz]
for arr in arrowlist:
ax.add_artist(arr)
# Rotate the view
ax.view_init(30, 30)
# Annotate the axes, shifts are ad-hoc for this (30, 30) view
xp, yp, _ = proj3d.proj_transform(arlen, 0, 0, ax.get_proj())
plt.annotate("x", xy=(xp, yp), xytext=(-3, -8), textcoords='offset points',
ha='right', va='bottom')
xp, yp, _ = proj3d.proj_transform(0, arlen, 0, ax.get_proj())
plt.annotate("y", xy=(xp, yp), xytext=(6, -5), textcoords='offset points',
ha='right', va='bottom')
xp, yp, _ = proj3d.proj_transform(0, 0, arlen, ax.get_proj())
plt.annotate("z", xy=(xp, yp), xytext=(2, 0), textcoords='offset points',
ha='right', va='bottom')
plt.title(title)
plt.show()
def plot_state_city(rho, title=""):
"""Plot the cityscape of quantum state.
Plot two 3d bargraphs (two dimenstional) of the mixed state rho
Args:
rho (np.array[[complex]]): array of dimensions 2**n x 2**nn complex
numbers
title (str): a string that represents the plot title
Returns:
none: plot is shown with matplotlib to the screen
"""
num = int(np.log2(len(rho)))
if basis is None:
basis = ['X', 'Y', 'Z']
elif isinstance(basis, dict):
basis = basis.keys()
# get the real and imag parts of rho
datareal = np.real(rho)
dataimag = np.imag(rho)
# get the labels
column_names = [bin(i)[2:].zfill(num) for i in range(2**num)]
row_names = [bin(i)[2:].zfill(num) for i in range(2**num)]
lx = len(datareal[0]) # Work out matrix dimensions
ly = len(datareal[:, 0])
xpos = np.arange(0, lx, 1) # Set up a mesh of positions
ypos = np.arange(0, ly, 1)
xpos, ypos = np.meshgrid(xpos+0.25, ypos+0.25)
xpos = xpos.flatten()
ypos = ypos.flatten()
zpos = np.zeros(lx*ly)
dx = 0.5 * np.ones_like(zpos) # width of bars
dy = dx.copy()
dzr = datareal.flatten()
dzi = dataimag.flatten()
fig = plt.figure(figsize=(8, 8))
ax1 = fig.add_subplot(2, 1, 1, projection='3d')
ax1.bar3d(xpos, ypos, zpos, dx, dy, dzr, color="g", alpha=0.5)
ax2 = fig.add_subplot(2, 1, 2, projection='3d')
ax2.bar3d(xpos, ypos, zpos, dx, dy, dzi, color="g", alpha=0.5)
ax1.set_xticks(np.arange(0.5, lx+0.5, 1))
ax1.set_yticks(np.arange(0.5, ly+0.5, 1))
ax1.axes.set_zlim3d(-1.0, 1.0001)
ax1.set_zticks(np.arange(-1, 1, 0.5))
ax1.w_xaxis.set_ticklabels(row_names, fontsize=12, rotation=45)
ax1.w_yaxis.set_ticklabels(column_names, fontsize=12, rotation=-22.5)
# ax1.set_xlabel('basis state', fontsize=12)
# ax1.set_ylabel('basis state', fontsize=12)
ax1.set_zlabel("Real[rho]")
ax2.set_xticks(np.arange(0.5, lx+0.5, 1))
ax2.set_yticks(np.arange(0.5, ly+0.5, 1))
ax2.axes.set_zlim3d(-1.0, 1.0001)
ax2.set_zticks(np.arange(-1, 1, 0.5))
ax2.w_xaxis.set_ticklabels(row_names, fontsize=12, rotation=45)
ax2.w_yaxis.set_ticklabels(column_names, fontsize=12, rotation=-22.5)
# ax2.set_xlabel('basis state', fontsize=12)
# ax2.set_ylabel('basis state', fontsize=12)
ax2.set_zlabel("Imag[rho]")
plt.title(title)
plt.show()
def plot_state_paulivec(rho, title=""):
"""Plot the paulivec representation of a quantum state.
Plot a bargraph of the mixed state rho over the pauli matricies
Args:
rho (np.array[[complex]]): array of dimensions 2**n x 2**nn complex
numbers
title (str): a string that represents the plot title
Returns:
none: plot is shown with matplotlib to the screen
"""
num = int(np.log2(len(rho)))
labels = list(map(lambda x: x.to_label(), pauli_group(num)))
values = list(map(lambda x: np.real(np.trace(np.dot(x.to_matrix(), rho))),
pauli_group(num)))
numelem = len(values)
ind = np.arange(numelem) # the x locations for the groups
width = 0.5 # the width of the bars
fig, ax = plt.subplots()
ax.grid(zorder=0)
ax.bar(ind, values, width, color='seagreen')
# add some text for labels, title, and axes ticks
ax.set_ylabel('Expectation value', fontsize=12)
ax.set_xticks(ind)
ax.set_yticks([-1, -0.5, 0, 0.5, 1])
ax.set_xticklabels(labels, fontsize=12, rotation=70)
ax.set_xlabel('Pauli', fontsize=12)
ax.set_ylim([-1, 1])
plt.title(title)
plt.show()
def n_choose_k(n, k):
"""Return the number of combinations for n choose k.
Args:
n (int): the total number of options .
k (int): The number of elements.
Returns:
int: returns the binomial coefficient
"""
if n == 0:
return 0
if not tomos:
return [{qubit: b} for b in basis]
else:
return reduce(lambda x, y: x * y[0] / y[1],
zip(range(n - k + 1, n + 1),
range(1, k + 1)), 1)
def lex_index(n, k, lst):
"""Return the lex index of a combination..
Args:
n (int): the total number of options .
k (int): The number of elements.
lst
Returns:
int: returns int index for lex order
ret = []
for t in tomos:
for b in basis:
t[qubit] = b
ret.append(t.copy())
return ret
"""
assert len(lst) == k, "list should have length k"
comb = list(map(lambda x: n - 1 - x, lst))
dualm = sum([n_choose_k(comb[k - 1 - i], i + 1) for i in range(k)])
return int(dualm)
if isinstance(qubits, int):
qubits = [qubits]
if basis is None:
basis = ['X', 'Y', 'Z']
elif isinstance(basis, dict):
basis = basis.keys()
return [dict(zip(qubits, b)) for b in product(basis, repeat=len(qubits))]
def bit_string_index(s):
"""Return the index of a string of 0s and 1s."""
n = len(s)
k = s.count("1")
assert s.count("0") == n - k, "s must be a string of 0 and 1"
ones = [pos for pos, char in enumerate(s) if char == "1"]
return lex_index(n, k, ones)
def phase_to_color_wheel(complex_number):
"""Map a phase of a complexnumber to a color in (r,g,b).
complex_number is phase is first mapped to angle in the range
[0, 2pi] and then to a color wheel with blue at zero phase.
def __add_meas_circuits(Q_program, name, qubits, qreg, creg, prep=None):
"""Helper function.
"""
angles = np.angle(complex_number)
angle_round = int(((angles + 2 * np.pi) % (2 * np.pi))/np.pi*6)
# print(angleround)
color_map = {0: (0, 0, 1), # blue,
1: (0.5, 0, 1), # blue-violet
2: (1, 0, 1), # violet
3: (1, 0, 0.5), # red-violet,
4: (1, 0, 0), # red
5: (1, 0.5, 0), # red-oranage,
6: (1, 1, 0), # orange
7: (0.5, 1, 0), # orange-yellow
8: (0, 1, 0), # yellow,
9: (0, 1, 0.5), # yellow-green,
10: (0, 1, 1), # green,
11: (0, 0.5, 1) # green-blue,
}
return color_map[angle_round]
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 = []
def plot_state_qsphere(rho):
"""Plot the qsphere representation of a quantum state."""
num = int(np.log2(len(rho)))
# get the eigenvectors and egivenvalues
we, stateall = la.eigh(rho)
for i in range(2**num):
# start with the max
probmix = we.max()
prob_location = we.argmax()
if probmix > 0.001:
print("The " + str(i) + "th eigenvalue = " + str(probmix))
# get the max eigenvalue
state = stateall[:, prob_location]
loc = np.absolute(state).argmax()
# get the element location closes to lowest bin representation.
for j in range(2**num):
test = np.absolute(np.absolute(state[j])
- np.absolute(state[loc]))
if test < 0.001:
loc = j
break
# remove the global phase
angles = (np.angle(state[loc]) + 2 * np.pi) % (2 * np.pi)
angleset = np.exp(-1j*angles)
# print(state)
# print(angles)
state = angleset*state
# print(state)
state.flatten()
# start the plotting
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.axes.set_xlim3d(-1.0, 1.0)
ax.axes.set_ylim3d(-1.0, 1.0)
ax.axes.set_zlim3d(-1.0, 1.0)
ax.set_aspect("equal")
ax.axes.grid(False)
# Plot semi-transparent sphere
u = np.linspace(0, 2 * np.pi, 25)
v = np.linspace(0, np.pi, 25)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, rstride=1, cstride=1, color='k',
alpha=0.05, linewidth=0)
# wireframe
# Get rid of the panes
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
for dic in __tomo_dicts(qubits):
# Get rid of the spines
ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
# Get rid of the ticks
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
# Construct meas circuit name
label = '_meas'
for qubit, op in dic.items():
label += op + str(qubit)
circ = Q_program.create_circuit(label, [qreg], [creg])
d = num
for i in range(2**num):
# get x,y,z points
element = bin(i)[2:].zfill(num)
weight = element.count("1")
zvalue = -2 * weight / d + 1
number_of_divisions = n_choose_k(d, weight)
weight_order = bit_string_index(element)
# if weight_order >= number_of_divisions / 2:
# com_key = compliment(element)
# weight_order_temp = bit_string_index(com_key)
# weight_order = np.floor(
# number_of_divisions / 2) + weight_order_temp + 1
angle = (weight_order) * 2 * np.pi / number_of_divisions
xvalue = np.sqrt(1 - zvalue**2) * np.cos(angle)
yvalue = np.sqrt(1 - zvalue**2) * np.sin(angle)
ax.plot([xvalue], [yvalue], [zvalue],
markerfacecolor=(.5, .5, .5),
markeredgecolor=(.5, .5, .5),
marker='o', markersize=10, alpha=1)
# get prob and angle - prob will be shade and angle color
prob = np.real(np.dot(state[i], state[i].conj()))
colorstate = phase_to_color_wheel(state[i])
a = Arrow3D([0, xvalue], [0, yvalue], [0, zvalue],
mutation_scale=20, alpha=prob, arrowstyle="-",
color=colorstate, lw=10)
ax.add_artist(a)
# add weight lines
for weight in range(d + 1):
theta = np.linspace(-2 * np.pi, 2 * np.pi, 100)
z = -2 * weight / d + 1
r = np.sqrt(1 - z**2)
x = r * np.cos(theta)
y = r * np.sin(theta)
ax.plot(x, y, z, color=(.5, .5, .5))
# add center point
ax.plot([0], [0], [0], markerfacecolor=(.5, .5, .5),
markeredgecolor=(.5, .5, .5), marker='o', markersize=10,
alpha=1)
plt.show()
we[prob_location] = 0
else:
break
# add gates to circuit
for qubit, op in dic.items():
circ.barrier(qreg[qubit])
if op == "X":
circ.h(qreg[qubit])
elif op == "Y":
circ.sdg(qreg[qubit])
circ.h(qreg[qubit])
if prep:
circ.barrier(qreg[qubit])
else:
circ.measure(qreg[qubit], creg[qubit])
# add circuit to QuantumProgram
Q_program.add_circuit(name+label, orig + circ)
# add label to output
labels.append(name+label)
# delete temp circuit
del Q_program._QuantumProgram__quantum_program['circuits'][label]
return labels
def plot_state(rho, method='city'):
"""Plot the quantum state."""
num = int(np.log2(len(rho)))
# Need updating to check its a matrix
if method == 'city':
plot_state_city(rho)
elif method == "paulivec":
plot_state_paulivec(rho)
elif method == "qsphere":
plot_state_qsphere(rho)
elif method == "bloch":
for i in range(num):
bloch_state = list(map(lambda x: np.real(np.trace(
np.dot(x.to_matrix(), rho))),
pauli_singles(i, num)))
plot_bloch_vector(bloch_state, "qubit " + str(i))
else:
print("No method given")
def __add_prep_circuits(Q_program, name, qubits, qreg, creg):
"""Helper function.
"""
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 = []
for dic in __tomo_dicts(qubits):
# Construct meas circuit name
label_p = '_prep' # prepare in positive eigenstate
label_m = '_prep' # prepare in negative eigenstate
for qubit, op in dic.items():
label_p += op + 'p' + str(qubit)
label_m+= op + 'm' + str(qubit)
circ_p = Q_program.create_circuit(label_p, [qreg], [creg])
circ_m = Q_program.create_circuit(label_m, [qreg], [creg])
# add gates to circuit
for qubit, op in dic.items():
if op == "X":
circ_p.h(qreg[qubit])
circ_m.x(qreg[qubit])
circ_m.h(qreg[qubit])
elif op == "Y":
circ_p.h(qreg[qubit])
circ_p.s(qreg[qubit])
circ_m.x(qreg[qubit])
circ_m.h(qreg[qubit])
circ_m.s(qreg[qubit])
elif op == "Z":
circ_m.x(qreg[qubit])
# add circuit to QuantumProgram
Q_program.add_circuit(name+label_p, circ_p + orig)
Q_program.add_circuit(name+label_m, circ_m + orig)
# add label to output
labels += [name+label_p, name+label_m]
# delete temp circuit
del Q_program._QuantumProgram__quantum_program['circuits'][label_p]
del Q_program._QuantumProgram__quantum_program['circuits'][label_m]
return labels
###############################################################
# Constructing tomographic measurement circuits
# Tomography circuit labels
###############################################################
def __tomo_labels(name, qubits, basis=None, subscript=None):
"""Helper function.
"""
if subscript is None:
subscript = ''
labels = []
for dic in __tomo_dicts(qubits, basis):
label = ''
for qubit, op in dic.items():
label += op + subscript + str(qubit)
labels.append(name+label)
return labels
def build_keys_helper(keys, qubit):
"""Return array of measurement strings ['Xj', 'Yj', 'Zj'] for qubit=j."""
tmp = []
for k in keys:
for b in ["X", "Y", "Z"]:
tmp.append(k + b + str(qubit))
return tmp
def state_tomography_labels(name, qubits, basis=None):
"""
"""
return __tomo_labels(name + '_meas', qubits, basis)
# Make private for now, since fit method not yet implemented
def __process_tomography_labels(name, qubits, basis=None):
"""
"""
preps = __tomo_labels(name + '_prep', qubits, basis, subscript='p')
preps += __tomo_labels(name + '_prep', qubits, basis, subscript='m')
return reduce(lambda acc, c: acc + __tomo_labels(c + '_meas', qubits, basis), preps, [])
###############################################################
# Tomography measurement outcomes
###############################################################
# 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 build_tomo_keys(circuit, qubit_list):
def __counts_keys(n):
"""helper function.
"""
For input circuit string returns an array of all measurement circits orded in
lexicographic order from last qubit to first qubit.
Example:
qubit_list = [0]: [circuitX0, circuitY0, circuitZ0].
qubit_list = [0,1]: [circuitX1X0, circuitX1Y0, circuitX1Z0, circuitY1X0,..., circuitZ1Z0].
"""
keys = [circuit]
for j in sorted(qubit_list, reverse=True):
keys = build_keys_helper(keys, j)
return keys
return [bin(j)[2:].zfill(n) for j in range(2 ** n)]
# We need to build circuits in QASM lexical order, not standard!
def build_tomo_circuit_helper(Q_program, circuits, qreg: str, creg: str, qubit: int):
"""
Adds measurements for the qubit=j to the input circuits, so if circuits = [c0, c1,...]
circuits-> [c0Xj, c0Yj, c0Zj, c1Xj,...]
"""
for c in circuits:
circ = Q_program.get_circuit(c)
for b in ["X","Y","Z"]:
meas = b+str(qubit)
tmp = Q_program.create_circuit(meas, [qreg],[creg])
qr = Q_program.get_quantum_registers(qreg)
cr = Q_program.get_classical_registers(creg)
if b == "X":
tmp.u2(0., np.pi, qr[qubit])
if b == "Y":
tmp.u2(0., np.pi / 2., qr[qubit])
tmp.measure(qr[qubit], cr[qubit])
Q_program.add_circuit(c+meas, circ + tmp)
def __get_basis_ops(tup, basis):
return reduce(lambda acc, b: [np.kron(a,j) for a in acc for j in basis[b]], tup, [1])
def build_tomo_circuits(Q_program, circuit, qreg: str, creg: str, qubit_list):
def __counts_basis(n, basis):
return [dict(zip(__counts_keys(n), __get_basis_ops(key, basis)))
for key in product(basis.keys(), repeat=n)]
def marginal_counts(counts, meas_qubits):
"""
Generates circuits in the input QuantumProgram for implementing complete quantum
state tomography of the state of the qubits in qubit_list prepared by circuit.
Returns a list of the measurement outcome strings for meas_qubits qubits in an
nq qubit system.
"""
circ = [circuit]
for j in sorted(qubit_list, reverse=True):
build_tomo_circuit_helper(Q_program, circ, qreg, creg, j)
circ = build_keys_helper(circ, j)
# Extract total number of qubits from count keys
nq = len(list(counts.keys())[0])
# keys for measured qubits only
qs = sorted(meas_qubits, reverse=True)
meas_keys = __counts_keys(len(qs))
# get regex match strings for suming outcomes of other qubits
rgx = [reduce(lambda x, y: (key[qs.index(y)] if y in qs else '\d')+x, range(nq), '') for key in meas_keys]
# build the return list
meas_counts = []
for m in rgx:
c = 0
for key, val in counts.items():
if match(m, key):
c += val
meas_counts.append(c)
# return as counts dict on measured qubits only
return dict(zip(meas_keys, meas_counts))
def state_tomography_data(Q_program, name, meas_qubits, backend=None, basis=None):
"""
"""
if basis is None:
basis = __default_basis
labels = state_tomography_labels(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]
meas_basis = __counts_basis(len(meas_qubits), basis)
ret = [{'counts': i, 'meas_basis': j, 'shots': k} for i, j, k in zip(counts, meas_basis, shots)]
return ret
###############################################################
@ -486,18 +298,7 @@ def build_tomo_circuits(Q_program, circuit, qreg: str, creg: str, qubit_list):
###############################################################
def vectorize(op):
"""Flatten an operator to a column-major vector."""
return op.flatten(order='F')
def devectorize(v):
"""Devectorize a column-major vectorized square matrix."""
d = int(np.sqrt(v.size))
return v.reshape(d, d, order='F')
def meas_basis_matrix(meas_basis):
def __tomo_basis_matrix(meas_basis):
"""Return a matrix of vectorized measurement operators.
Measurement operators S = sum_j |j><M_j|.
@ -507,148 +308,91 @@ def meas_basis_matrix(meas_basis):
S = np.array([vectorize(m).conj() for m in meas_basis])
return S.reshape(n, d)
def outer(v1, v2=None):
"""
Returns the matrix |v1><v2| resulting from the outer product of two vectors.
"""
if v2 is None:
u = v1.conj()
else:
u = v2.conj()
return np.outer(v1, u)
def wizard(rho, normalize_flag = True, epsilon = 0.):
def __tomo_linear_inv(freqs, ops, weights=None, trace=1):
"""
"""
# get weights matrix
if weights is not None:
W = np.array(weights)
if W.ndim == 1:
W = np.diag(W)
# Get basis S matrix
S = np.array([vectorize(m).conj() for m in ops]).reshape(len(ops), ops[0].size)
if weights is not None:
S = np.dot(W, S) # W.S
# get frequencies vec
v = np.array(freqs) # |n>
if weights is not None:
v = np.dot(W, freqs) #W.|n>
Sdg = S.T.conj() # S^*.W^*
inv = np.linalg.pinv(np.dot(Sdg,S))
# linear inversion of freqs
ret = devectorize(np.dot(inv, np.dot(Sdg, v)))
# renormalize to input trace value
ret = trace * ret / np.trace(ret)
return ret
def __state_leastsq_fit(state_data, weights=None, beta=0.5):
"""
"""
ks = state_data[0]['counts'].keys()
# Get counts and shots
ns = np.array([dat['counts'][k] for dat in state_data for k in ks])
shots = np.array([dat['shots'] for dat in state_data for k in ks])
# convert to freqs
freqs = (ns + beta) / (shots + 2 * beta)
# Use standard least squares fitting weights
if weights is None:
weights = np.sqrt(shots / (freqs * (1 - freqs)))
# Get measurement basis ops
ops = [dat['meas_basis'][k] for dat in state_data for k in ks]
return __tomo_linear_inv(freqs, ops, weights, trace=1)
def __wizard(rho, epsilon=0.):
"""
Maps an operator to the nearest positive semidefinite operator
by setting negative eigenvalues to zero and rescaling the positive
eigenvalues.
See arXiv:1106.5458 [quant-ph]
"""
#print("Using wizard method to constrain positivity")
if normalize_flag:
rho = rho / np.trace(rho)
dim = len(rho)
rho_wizard = np.zeros([dim, dim])
v, w = np.linalg.eigh(rho) # v eigenvecrors v[0] < v[1] <...
v, w = np.linalg.eigh(rho) # v eigenvecrors v[0] < v[1] <...
for j in range(dim):
if v[j] < epsilon:
tmp = v[j]
v[j] = 0.
# redistribute loop
x = 0.
for k in range(j+1,dim):
for k in range(j + 1, dim):
x += tmp / (dim-(j+1))
v[k] = v[k] + tmp / (dim -(j+1))
v[k] = v[k] + tmp / (dim - (j+1))
for j in range(dim):
rho_wizard = rho_wizard + v[j] * outer(w[:,j])
rho_wizard = rho_wizard + v[j] * outer(w[:, j])
return rho_wizard
def fit_state(freqs, meas_basis, weights=None, normalize_flag = True, wizard_flag = False):
def fit_state(state_data, method='wizard', options=None):
"""
Returns a matrix reconstructed by unconstrained least-squares fitting.
"""
if weights is None:
W = np.eye(len(freqs)) # use uniform weights
if method in ['wizard', 'leastsq']:
rho = __state_leastsq_fit(state_data)
if method == 'wizard':
rho = __wizard(rho)
else:
W = np.array(np.diag(weights))
S = np.dot(W, meas_basis_matrix(meas_basis)) # actually W.S
v = np.dot(W, freqs) # W|f>
v = np.array(np.dot(S.T.conj(), v)) # S^*.W^*W.|f>
inv = np.linalg.pinv(np.dot(S.T.conj(), S)) # pseudo inverse of S^*.W^*.W.S
v = np.dot(inv, v) # |rho>
rho = devectorize(v)
if normalize_flag:
rho = rho / np.trace(rho)
if wizard_flag:
#rho_wizard = wizard(rho,normalize_flag)
rho = wizard(rho, normalize_flag=normalize_flag)
return rho
###############################################################
# Parsing measurement data for reconstruction
###############################################################
def nqubit_basis(n):
"""
Returns the measurement basis for n-qubits in the correct order for the
meas_outcome_strings function.
"""
b1 = np.array([
np.array([[0.5, 0.5], [0.5, 0.5]]), np.array([[0.5, -0.5], [-0.5, 0.5]]), # Xp, Xm
np.array([[0.5, -0.5j], [0.5j, 0.5]]), np.array([[0.5, 0.5j], [-0.5j, 0.5]]), # Yp, Ym
np.array([[1, 0], [0, 0]]), np.array([[0, 0], [0, 1]]), # Zp, Zm
])
if n == 1:
return b1
else:
bnm1 = nqubit_basis(n-1)
m = 2**(n-1)
d = bnm1.size // m**3
return np.kron( bnm1.reshape(d, m, m, m), b1.reshape(3, 2, 2, 2)).reshape(3*d * 2*m, 2*m, 2*m)
def meas_outcome_strings(nq):
"""
Returns a list of the measurement outcome strings for nq qubits.
"""
return [bin(j)[2:].zfill(nq) for j in range(2**nq)]
def tomo_outcome_strings(meas_qubits,nq=None):
"""
Returns a list of the measurement outcome strings for meas_qubits qubits in an
nq qubit system.
"""
if nq is None:
nq = len(meas_qubits)
qs = sorted(meas_qubits, reverse=True)
outs = meas_outcome_strings(len(qs))
res = [];
for s in outs:
label = ""
for j in range(nq):
if j in qs:
label = s[qs.index(j)] + label
else:
label = str(0) + label
res.append(label)
return res
def none_to_zero(val):
"""
Returns 0 if the input argument is None, else it returns the input.
"""
if val is None:
return 0
else:
return val
###############################################################
# Putting it all together
###############################################################
def state_tomography(Q_program, tomo_circuits, shots, nq,
meas_qubits, method='leastsq_wizard'):
"""
Returns the reconstructed density matrix.
"""
m = len(meas_qubits)
counts = np.array([none_to_zero(Q_program.get_counts(c).get(s))
for c in tomo_circuits
for s in tomo_outcome_strings(meas_qubits, nq)])
if method == 'leastsq_wizard':
return fit_state(counts / shots, nqubit_basis(m), normalize_flag=True, wizard_flag=True)
elif method == 'least_sq':
return fit_state(counts / shots, nqubit_basis(m), normalize_flag=True, wizard_flag=False)
else:
print("error: unknown reconstruction method")
def state_fidelity(rho, psi):
"""Return the state fidelity.
F between a density matrix rho and a target pure state psi.
"""
return np.sqrt(np.abs(np.dot(psi, np.dot(rho, psi))))
def is_pos_def(x):
return np.all(np.linalg.eigvals(x) > 0)
# TODO: raise exception for unknown method
print('error: method unrecongnized')
pass
return rho

413
tools/vizualization.py Normal file
View File

@ -0,0 +1,413 @@
# -*- coding: utf-8 -*-
# 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
# =============================================================================
"""
Visualization functions for quantum states.
Author: Jay Gambetta
Andrew Cross
Christopher J. Wood
"""
import numpy as np
from functools import reduce
from scipy import linalg as la
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch
from tools.pauli import pauli_group, pauli_singles
###############################################################
# Plotting tools
###############################################################
class Arrow3D(FancyArrowPatch):
"""Standard 3D arrow."""
def __init__(self, xs, ys, zs, *args, **kwargs):
"""Create arrow."""
FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
self._verts3d = xs, ys, zs
def draw(self, renderer):
"""Draw the arrow."""
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
FancyArrowPatch.draw(self, renderer)
def plot_bloch_vector(bloch, title=""):
"""Plot the Bloch sphere.
Plot a sphere, axes, the Bloch vector, and its projections onto each axis.
Args:
bloch (list[double]): array of three elements where [<x>, <y>,<z>]
title (str): a string that represents the plot title
Returns:
none: plot is shown with matplotlib to the screen
"""
# Set arrow lengths
arlen = 1.3
# Plot semi-transparent sphere
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='3d')
ax.set_aspect("equal")
ax.plot_surface(x, y, z, color=(.5, .5, .5), alpha=0.1)
# Plot arrows (axes, Bloch vector, its projections)
xa = Arrow3D([0, arlen], [0, 0], [0, 0], mutation_scale=20, lw=1,
arrowstyle="-|>", color=(.5, .5, .5))
ya = Arrow3D([0, 0], [0, arlen], [0, 0], mutation_scale=20, lw=1,
arrowstyle="-|>", color=(.5, .5, .5))
za = Arrow3D([0, 0], [0, 0], [0, arlen], mutation_scale=20, lw=1,
arrowstyle="-|>", color=(.5, .5, .5))
a = Arrow3D([0, bloch[0]], [0, bloch[1]], [0, bloch[2]], mutation_scale=20,
lw=2, arrowstyle="simple", color="k")
bax = Arrow3D([0, bloch[0]], [0, 0], [0, 0], mutation_scale=20, lw=2,
arrowstyle="-", color="r")
bay = Arrow3D([0, 0], [0, bloch[1]], [0, 0], mutation_scale=20, lw=2,
arrowstyle="-", color="g")
baz = Arrow3D([0, 0], [0, 0], [0, bloch[2]], mutation_scale=20, lw=2,
arrowstyle="-", color="b")
arrowlist = [xa, ya, za, a, bax, bay, baz]
for arr in arrowlist:
ax.add_artist(arr)
# Rotate the view
ax.view_init(30, 30)
# Annotate the axes, shifts are ad-hoc for this (30, 30) view
xp, yp, _ = proj3d.proj_transform(arlen, 0, 0, ax.get_proj())
plt.annotate("x", xy=(xp, yp), xytext=(-3, -8), textcoords='offset points',
ha='right', va='bottom')
xp, yp, _ = proj3d.proj_transform(0, arlen, 0, ax.get_proj())
plt.annotate("y", xy=(xp, yp), xytext=(6, -5), textcoords='offset points',
ha='right', va='bottom')
xp, yp, _ = proj3d.proj_transform(0, 0, arlen, ax.get_proj())
plt.annotate("z", xy=(xp, yp), xytext=(2, 0), textcoords='offset points',
ha='right', va='bottom')
plt.title(title)
plt.show()
def plot_state_city(rho, title=""):
"""Plot the cityscape of quantum state.
Plot two 3d bargraphs (two dimenstional) of the mixed state rho
Args:
rho (np.array[[complex]]): array of dimensions 2**n x 2**nn complex
numbers
title (str): a string that represents the plot title
Returns:
none: plot is shown with matplotlib to the screen
"""
num = int(np.log2(len(rho)))
# get the real and imag parts of rho
datareal = np.real(rho)
dataimag = np.imag(rho)
# get the labels
column_names = [bin(i)[2:].zfill(num) for i in range(2**num)]
row_names = [bin(i)[2:].zfill(num) for i in range(2**num)]
lx = len(datareal[0]) # Work out matrix dimensions
ly = len(datareal[:, 0])
xpos = np.arange(0, lx, 1) # Set up a mesh of positions
ypos = np.arange(0, ly, 1)
xpos, ypos = np.meshgrid(xpos+0.25, ypos+0.25)
xpos = xpos.flatten()
ypos = ypos.flatten()
zpos = np.zeros(lx*ly)
dx = 0.5 * np.ones_like(zpos) # width of bars
dy = dx.copy()
dzr = datareal.flatten()
dzi = dataimag.flatten()
fig = plt.figure(figsize=(8, 8))
ax1 = fig.add_subplot(2, 1, 1, projection='3d')
ax1.bar3d(xpos, ypos, zpos, dx, dy, dzr, color="g", alpha=0.5)
ax2 = fig.add_subplot(2, 1, 2, projection='3d')
ax2.bar3d(xpos, ypos, zpos, dx, dy, dzi, color="g", alpha=0.5)
ax1.set_xticks(np.arange(0.5, lx+0.5, 1))
ax1.set_yticks(np.arange(0.5, ly+0.5, 1))
ax1.axes.set_zlim3d(-1.0, 1.0001)
ax1.set_zticks(np.arange(-1, 1, 0.5))
ax1.w_xaxis.set_ticklabels(row_names, fontsize=12, rotation=45)
ax1.w_yaxis.set_ticklabels(column_names, fontsize=12, rotation=-22.5)
# ax1.set_xlabel('basis state', fontsize=12)
# ax1.set_ylabel('basis state', fontsize=12)
ax1.set_zlabel("Real[rho]")
ax2.set_xticks(np.arange(0.5, lx+0.5, 1))
ax2.set_yticks(np.arange(0.5, ly+0.5, 1))
ax2.axes.set_zlim3d(-1.0, 1.0001)
ax2.set_zticks(np.arange(-1, 1, 0.5))
ax2.w_xaxis.set_ticklabels(row_names, fontsize=12, rotation=45)
ax2.w_yaxis.set_ticklabels(column_names, fontsize=12, rotation=-22.5)
# ax2.set_xlabel('basis state', fontsize=12)
# ax2.set_ylabel('basis state', fontsize=12)
ax2.set_zlabel("Imag[rho]")
plt.title(title)
plt.show()
def plot_state_paulivec(rho, title=""):
"""Plot the paulivec representation of a quantum state.
Plot a bargraph of the mixed state rho over the pauli matricies
Args:
rho (np.array[[complex]]): array of dimensions 2**n x 2**nn complex
numbers
title (str): a string that represents the plot title
Returns:
none: plot is shown with matplotlib to the screen
"""
num = int(np.log2(len(rho)))
labels = list(map(lambda x: x.to_label(), pauli_group(num)))
values = list(map(lambda x: np.real(np.trace(np.dot(x.to_matrix(), rho))),
pauli_group(num)))
numelem = len(values)
ind = np.arange(numelem) # the x locations for the groups
width = 0.5 # the width of the bars
fig, ax = plt.subplots()
ax.grid(zorder=0)
ax.bar(ind, values, width, color='seagreen')
# add some text for labels, title, and axes ticks
ax.set_ylabel('Expectation value', fontsize=12)
ax.set_xticks(ind)
ax.set_yticks([-1, -0.5, 0, 0.5, 1])
ax.set_xticklabels(labels, fontsize=12, rotation=70)
ax.set_xlabel('Pauli', fontsize=12)
ax.set_ylim([-1, 1])
plt.title(title)
plt.show()
def n_choose_k(n, k):
"""Return the number of combinations for n choose k.
Args:
n (int): the total number of options .
k (int): The number of elements.
Returns:
int: returns the binomial coefficient
"""
if n == 0:
return 0
else:
return reduce(lambda x, y: x * y[0] / y[1],
zip(range(n - k + 1, n + 1),
range(1, k + 1)), 1)
def lex_index(n, k, lst):
"""Return the lex index of a combination..
Args:
n (int): the total number of options .
k (int): The number of elements.
lst
Returns:
int: returns int index for lex order
"""
assert len(lst) == k, "list should have length k"
comb = list(map(lambda x: n - 1 - x, lst))
dualm = sum([n_choose_k(comb[k - 1 - i], i + 1) for i in range(k)])
return int(dualm)
def bit_string_index(s):
"""Return the index of a string of 0s and 1s."""
n = len(s)
k = s.count("1")
assert s.count("0") == n - k, "s must be a string of 0 and 1"
ones = [pos for pos, char in enumerate(s) if char == "1"]
return lex_index(n, k, ones)
def phase_to_color_wheel(complex_number):
"""Map a phase of a complexnumber to a color in (r,g,b).
complex_number is phase is first mapped to angle in the range
[0, 2pi] and then to a color wheel with blue at zero phase.
"""
angles = np.angle(complex_number)
angle_round = int(((angles + 2 * np.pi) % (2 * np.pi))/np.pi*6)
# print(angleround)
color_map = {0: (0, 0, 1), # blue,
1: (0.5, 0, 1), # blue-violet
2: (1, 0, 1), # violet
3: (1, 0, 0.5), # red-violet,
4: (1, 0, 0), # red
5: (1, 0.5, 0), # red-oranage,
6: (1, 1, 0), # orange
7: (0.5, 1, 0), # orange-yellow
8: (0, 1, 0), # yellow,
9: (0, 1, 0.5), # yellow-green,
10: (0, 1, 1), # green,
11: (0, 0.5, 1) # green-blue,
}
return color_map[angle_round]
def plot_state_qsphere(rho):
"""Plot the qsphere representation of a quantum state."""
num = int(np.log2(len(rho)))
# get the eigenvectors and egivenvalues
we, stateall = la.eigh(rho)
for i in range(2**num):
# start with the max
probmix = we.max()
prob_location = we.argmax()
if probmix > 0.001:
print("The " + str(i) + "th eigenvalue = " + str(probmix))
# get the max eigenvalue
state = stateall[:, prob_location]
loc = np.absolute(state).argmax()
# get the element location closes to lowest bin representation.
for j in range(2**num):
test = np.absolute(np.absolute(state[j])
- np.absolute(state[loc]))
if test < 0.001:
loc = j
break
# remove the global phase
angles = (np.angle(state[loc]) + 2 * np.pi) % (2 * np.pi)
angleset = np.exp(-1j*angles)
# print(state)
# print(angles)
state = angleset*state
# print(state)
state.flatten()
# start the plotting
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.axes.set_xlim3d(-1.0, 1.0)
ax.axes.set_ylim3d(-1.0, 1.0)
ax.axes.set_zlim3d(-1.0, 1.0)
ax.set_aspect("equal")
ax.axes.grid(False)
# Plot semi-transparent sphere
u = np.linspace(0, 2 * np.pi, 25)
v = np.linspace(0, np.pi, 25)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, rstride=1, cstride=1, color='k',
alpha=0.05, linewidth=0)
# wireframe
# Get rid of the panes
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
# Get rid of the spines
ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
# Get rid of the ticks
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
d = num
for i in range(2**num):
# get x,y,z points
element = bin(i)[2:].zfill(num)
weight = element.count("1")
zvalue = -2 * weight / d + 1
number_of_divisions = n_choose_k(d, weight)
weight_order = bit_string_index(element)
# if weight_order >= number_of_divisions / 2:
# com_key = compliment(element)
# weight_order_temp = bit_string_index(com_key)
# weight_order = np.floor(
# number_of_divisions / 2) + weight_order_temp + 1
angle = (weight_order) * 2 * np.pi / number_of_divisions
xvalue = np.sqrt(1 - zvalue**2) * np.cos(angle)
yvalue = np.sqrt(1 - zvalue**2) * np.sin(angle)
ax.plot([xvalue], [yvalue], [zvalue],
markerfacecolor=(.5, .5, .5),
markeredgecolor=(.5, .5, .5),
marker='o', markersize=10, alpha=1)
# get prob and angle - prob will be shade and angle color
prob = np.real(np.dot(state[i], state[i].conj()))
colorstate = phase_to_color_wheel(state[i])
a = Arrow3D([0, xvalue], [0, yvalue], [0, zvalue],
mutation_scale=20, alpha=prob, arrowstyle="-",
color=colorstate, lw=10)
ax.add_artist(a)
# add weight lines
for weight in range(d + 1):
theta = np.linspace(-2 * np.pi, 2 * np.pi, 100)
z = -2 * weight / d + 1
r = np.sqrt(1 - z**2)
x = r * np.cos(theta)
y = r * np.sin(theta)
ax.plot(x, y, z, color=(.5, .5, .5))
# add center point
ax.plot([0], [0], [0], markerfacecolor=(.5, .5, .5),
markeredgecolor=(.5, .5, .5), marker='o', markersize=10,
alpha=1)
plt.show()
we[prob_location] = 0
else:
break
def plot_state(rho, method='city'):
"""Plot the quantum state."""
num = int(np.log2(len(rho)))
# Need updating to check its a matrix
if method == 'city':
plot_state_city(rho)
elif method == "paulivec":
plot_state_paulivec(rho)
elif method == "qsphere":
plot_state_qsphere(rho)
elif method == "bloch":
for i in range(num):
bloch_state = list(map(lambda x: np.real(np.trace(
np.dot(x.to_matrix(), rho))),
pauli_singles(i, num)))
plot_bloch_vector(bloch_state, "qubit " + str(i))
else:
print("No method given")

458
tutorial/sections/process_tomography_testing.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

598
tutorial/sections/state_tomography.ipynb
View File

File diff suppressed because one or more lines are too long