"""Qubits for quantum computing.
Todo:
* Finish implementing measurement logic. This should include POVM.
* Update docstrings.
* Update tests.
"""
from __future__ import print_function, division
import math
from sympy import Integer, log, Mul, Add, Pow, conjugate
from sympy.core.basic import sympify
from sympy.core.compatibility import string_types, range, SYMPY_INTS
from sympy.matrices import Matrix, zeros
from sympy.printing.pretty.stringpict import prettyForm
from sympy.physics.quantum.hilbert import ComplexSpace
from sympy.physics.quantum.state import Ket, Bra, State
from sympy.physics.quantum.qexpr import QuantumError
from sympy.physics.quantum.represent import represent
from sympy.physics.quantum.matrixutils import (
numpy_ndarray, scipy_sparse_matrix
)
from mpmath.libmp.libintmath import bitcount
__all__ = [
'Qubit',
'QubitBra',
'IntQubit',
'IntQubitBra',
'qubit_to_matrix',
'matrix_to_qubit',
'matrix_to_density',
'measure_all',
'measure_partial',
'measure_partial_oneshot',
'measure_all_oneshot'
]
#-----------------------------------------------------------------------------
# Qubit Classes
#-----------------------------------------------------------------------------
class QubitState(State):
"""Base class for Qubit and QubitBra."""
#-------------------------------------------------------------------------
# Initialization/creation
#-------------------------------------------------------------------------
@classmethod
def _eval_args(cls, args):
# If we are passed a QubitState or subclass, we just take its qubit
# values directly.
if len(args) == 1 and isinstance(args[0], QubitState):
return args[0].qubit_values
# Turn strings into tuple of strings
if len(args) == 1 and isinstance(args[0], string_types):
args = tuple(args[0])
args = sympify(args)
# Validate input (must have 0 or 1 input)
for element in args:
if not (element == 1 or element == 0):
raise ValueError(
"Qubit values must be 0 or 1, got: %r" % element)
return args
@classmethod
def _eval_hilbert_space(cls, args):
return ComplexSpace(2)**len(args)
#-------------------------------------------------------------------------
# Properties
#-------------------------------------------------------------------------
@property
def dimension(self):
"""The number of Qubits in the state."""
return len(self.qubit_values)
@property
def nqubits(self):
return self.dimension
@property
def qubit_values(self):
"""Returns the values of the qubits as a tuple."""
return self.label
#-------------------------------------------------------------------------
# Special methods
#-------------------------------------------------------------------------
def __len__(self):
return self.dimension
def __getitem__(self, bit):
return self.qubit_values[int(self.dimension - bit - 1)]
#-------------------------------------------------------------------------
# Utility methods
#-------------------------------------------------------------------------
def flip(self, *bits):
"""Flip the bit(s) given."""
newargs = list(self.qubit_values)
for i in bits:
bit = int(self.dimension - i - 1)
if newargs[bit] == 1:
newargs[bit] = 0
else:
newargs[bit] = 1
return self.__class__(*tuple(newargs))
[docs]class Qubit(QubitState, Ket):
"""A multi-qubit ket in the computational (z) basis.
We use the normal convention that the least significant qubit is on the
right, so ``|00001>`` has a 1 in the least significant qubit.
Parameters
==========
values : list, str
The qubit values as a list of ints ([0,0,0,1,1,]) or a string ('011').
Examples
========
Create a qubit in a couple of different ways and look at their attributes:
>>> from sympy.physics.quantum.qubit import Qubit
>>> Qubit(0,0,0)
|000>
>>> q = Qubit('0101')
>>> q
|0101>
>>> q.nqubits
4
>>> len(q)
4
>>> q.dimension
4
>>> q.qubit_values
(0, 1, 0, 1)
We can flip the value of an individual qubit:
>>> q.flip(1)
|0111>
We can take the dagger of a Qubit to get a bra:
>>> from sympy.physics.quantum.dagger import Dagger
>>> Dagger(q)
<0101|
>>> type(Dagger(q))
<class 'sympy.physics.quantum.qubit.QubitBra'>
Inner products work as expected:
>>> ip = Dagger(q)*q
>>> ip
<0101|0101>
>>> ip.doit()
1
"""
@classmethod
def dual_class(self):
return QubitBra
def _eval_innerproduct_QubitBra(self, bra, **hints):
if self.label == bra.label:
return Integer(1)
else:
return Integer(0)
def _represent_default_basis(self, **options):
return self._represent_ZGate(None, **options)
def _represent_ZGate(self, basis, **options):
"""Represent this qubits in the computational basis (ZGate).
"""
format = options.get('format', 'sympy')
n = 1
definite_state = 0
for it in reversed(self.qubit_values):
definite_state += n*it
n = n*2
result = [0]*(2**self.dimension)
result[int(definite_state)] = 1
if format == 'sympy':
return Matrix(result)
elif format == 'numpy':
import numpy as np
return np.matrix(result, dtype='complex').transpose()
elif format == 'scipy.sparse':
from scipy import sparse
return sparse.csr_matrix(result, dtype='complex').transpose()
def _eval_trace(self, bra, **kwargs):
indices = kwargs.get('indices', [])
#sort index list to begin trace from most-significant
#qubit
sorted_idx = list(indices)
if len(sorted_idx) == 0:
sorted_idx = list(range(0, self.nqubits))
sorted_idx.sort()
#trace out for each of index
new_mat = self*bra
for i in range(len(sorted_idx) - 1, -1, -1):
# start from tracing out from leftmost qubit
new_mat = self._reduced_density(new_mat, int(sorted_idx[i]))
if (len(sorted_idx) == self.nqubits):
#in case full trace was requested
return new_mat[0]
else:
return matrix_to_density(new_mat)
def _reduced_density(self, matrix, qubit, **options):
"""Compute the reduced density matrix by tracing out one qubit.
The qubit argument should be of type python int, since it is used
in bit operations
"""
def find_index_that_is_projected(j, k, qubit):
bit_mask = 2**qubit - 1
return ((j >> qubit) << (1 + qubit)) + (j & bit_mask) + (k << qubit)
old_matrix = represent(matrix, **options)
old_size = old_matrix.cols
#we expect the old_size to be even
new_size = old_size//2
new_matrix = Matrix().zeros(new_size)
for i in range(new_size):
for j in range(new_size):
for k in range(2):
col = find_index_that_is_projected(j, k, qubit)
row = find_index_that_is_projected(i, k, qubit)
new_matrix[i, j] += old_matrix[row, col]
return new_matrix
[docs]class QubitBra(QubitState, Bra):
"""A multi-qubit bra in the computational (z) basis.
We use the normal convention that the least significant qubit is on the
right, so ``|00001>`` has a 1 in the least significant qubit.
Parameters
==========
values : list, str
The qubit values as a list of ints ([0,0,0,1,1,]) or a string ('011').
See also
========
Qubit: Examples using qubits
"""
@classmethod
def dual_class(self):
return Qubit
class IntQubitState(QubitState):
"""A base class for qubits that work with binary representations."""
@classmethod
def _eval_args(cls, args, nqubits=None):
# The case of a QubitState instance
if len(args) == 1 and isinstance(args[0], QubitState):
return QubitState._eval_args(args)
# otherwise, args should be integer
elif not all((isinstance(a, (int, Integer)) for a in args)):
raise ValueError('values must be integers, got (%s)' % (tuple(type(a) for a in args),))
# use nqubits if specified
if nqubits is not None:
if not isinstance(nqubits, (int, Integer)):
raise ValueError('nqubits must be an integer, got (%s)' % type(nqubits))
if len(args) != 1:
raise ValueError(
'too many positional arguments (%s). should be (number, nqubits=n)' % (args,))
return cls._eval_args_with_nqubits(args[0], nqubits)
# For a single argument, we construct the binary representation of
# that integer with the minimal number of bits.
if len(args) == 1 and args[0] > 1:
#rvalues is the minimum number of bits needed to express the number
rvalues = reversed(range(bitcount(abs(args[0]))))
qubit_values = [(args[0] >> i) & 1 for i in rvalues]
return QubitState._eval_args(qubit_values)
# For two numbers, the second number is the number of bits
# on which it is expressed, so IntQubit(0,5) == |00000>.
elif len(args) == 2 and args[1] > 1:
return cls._eval_args_with_nqubits(args[0], args[1])
else:
return QubitState._eval_args(args)
@classmethod
def _eval_args_with_nqubits(cls, number, nqubits):
need = bitcount(abs(number))
if nqubits < need:
raise ValueError(
'cannot represent %s with %s bits' % (number, nqubits))
qubit_values = [(number >> i) & 1 for i in reversed(range(nqubits))]
return QubitState._eval_args(qubit_values)
def as_int(self):
"""Return the numerical value of the qubit."""
number = 0
n = 1
for i in reversed(self.qubit_values):
number += n*i
n = n << 1
return number
def _print_label(self, printer, *args):
return str(self.as_int())
def _print_label_pretty(self, printer, *args):
label = self._print_label(printer, *args)
return prettyForm(label)
_print_label_repr = _print_label
_print_label_latex = _print_label
[docs]class IntQubit(IntQubitState, Qubit):
"""A qubit ket that store integers as binary numbers in qubit values.
The differences between this class and ``Qubit`` are:
* The form of the constructor.
* The qubit values are printed as their corresponding integer, rather
than the raw qubit values. The internal storage format of the qubit
values in the same as ``Qubit``.
Parameters
==========
values : int, tuple
If a single argument, the integer we want to represent in the qubit
values. This integer will be represented using the fewest possible
number of qubits.
If a pair of integers and the second value is more than one, the first
integer gives the integer to represent in binary form and the second
integer gives the number of qubits to use.
List of zeros and ones is also accepted to generate qubit by bit pattern.
nqubits : int
The integer that represents the number of qubits.
This number should be passed with keyword ``nqubits=N``.
You can use this in order to avoid ambiguity of Qubit-style tuple of bits.
Please see the example below for more details.
Examples
========
Create a qubit for the integer 5:
>>> from sympy.physics.quantum.qubit import IntQubit
>>> from sympy.physics.quantum.qubit import Qubit
>>> q = IntQubit(5)
>>> q
|5>
We can also create an ``IntQubit`` by passing a ``Qubit`` instance.
>>> q = IntQubit(Qubit('101'))
>>> q
|5>
>>> q.as_int()
5
>>> q.nqubits
3
>>> q.qubit_values
(1, 0, 1)
We can go back to the regular qubit form.
>>> Qubit(q)
|101>
Please note that ``IntQubit`` also accepts a ``Qubit``-style list of bits.
So, the code below yields qubits 3, not a single bit ``1``.
>>> IntQubit(1, 1)
|3>
To avoid ambiguity, use ``nqubits`` parameter.
Use of this keyword is recommended especially when you provide the values by variables.
>>> IntQubit(1, nqubits=1)
|1>
>>> a = 1
>>> IntQubit(a, nqubits=1)
|1>
"""
@classmethod
def dual_class(self):
return IntQubitBra
def _eval_innerproduct_IntQubitBra(self, bra, **hints):
return Qubit._eval_innerproduct_QubitBra(self, bra)
[docs]class IntQubitBra(IntQubitState, QubitBra):
"""A qubit bra that store integers as binary numbers in qubit values."""
@classmethod
def dual_class(self):
return IntQubit
#-----------------------------------------------------------------------------
# Qubit <---> Matrix conversion functions
#-----------------------------------------------------------------------------
[docs]def matrix_to_qubit(matrix):
"""Convert from the matrix repr. to a sum of Qubit objects.
Parameters
----------
matrix : Matrix, numpy.matrix, scipy.sparse
The matrix to build the Qubit representation of. This works with
sympy matrices, numpy matrices and scipy.sparse sparse matrices.
Examples
========
Represent a state and then go back to its qubit form:
>>> from sympy.physics.quantum.qubit import matrix_to_qubit, Qubit
>>> from sympy.physics.quantum.gate import Z
>>> from sympy.physics.quantum.represent import represent
>>> q = Qubit('01')
>>> matrix_to_qubit(represent(q))
|01>
"""
# Determine the format based on the type of the input matrix
format = 'sympy'
if isinstance(matrix, numpy_ndarray):
format = 'numpy'
if isinstance(matrix, scipy_sparse_matrix):
format = 'scipy.sparse'
# Make sure it is of correct dimensions for a Qubit-matrix representation.
# This logic should work with sympy, numpy or scipy.sparse matrices.
if matrix.shape[0] == 1:
mlistlen = matrix.shape[1]
nqubits = log(mlistlen, 2)
ket = False
cls = QubitBra
elif matrix.shape[1] == 1:
mlistlen = matrix.shape[0]
nqubits = log(mlistlen, 2)
ket = True
cls = Qubit
else:
raise QuantumError(
'Matrix must be a row/column vector, got %r' % matrix
)
if not isinstance(nqubits, Integer):
raise QuantumError('Matrix must be a row/column vector of size '
'2**nqubits, got: %r' % matrix)
# Go through each item in matrix, if element is non-zero, make it into a
# Qubit item times the element.
result = 0
for i in range(mlistlen):
if ket:
element = matrix[i, 0]
else:
element = matrix[0, i]
if format == 'numpy' or format == 'scipy.sparse':
element = complex(element)
if element != 0.0:
# Form Qubit array; 0 in bit-locations where i is 0, 1 in
# bit-locations where i is 1
qubit_array = [int(i & (1 << x) != 0) for x in range(nqubits)]
qubit_array.reverse()
result = result + element*cls(*qubit_array)
# If sympy simplified by pulling out a constant coefficient, undo that.
if isinstance(result, (Mul, Add, Pow)):
result = result.expand()
return result
[docs]def matrix_to_density(mat):
"""
Works by finding the eigenvectors and eigenvalues of the matrix.
We know we can decompose rho by doing:
sum(EigenVal*|Eigenvect><Eigenvect|)
"""
from sympy.physics.quantum.density import Density
eigen = mat.eigenvects()
args = [[matrix_to_qubit(Matrix(
[vector, ])), x[0]] for x in eigen for vector in x[2] if x[0] != 0]
if (len(args) == 0):
return 0
else:
return Density(*args)
[docs]def qubit_to_matrix(qubit, format='sympy'):
"""Converts an Add/Mul of Qubit objects into it's matrix representation
This function is the inverse of ``matrix_to_qubit`` and is a shorthand
for ``represent(qubit)``.
"""
return represent(qubit, format=format)
#-----------------------------------------------------------------------------
# Measurement
#-----------------------------------------------------------------------------
[docs]def measure_all(qubit, format='sympy', normalize=True):
"""Perform an ensemble measurement of all qubits.
Parameters
==========
qubit : Qubit, Add
The qubit to measure. This can be any Qubit or a linear combination
of them.
format : str
The format of the intermediate matrices to use. Possible values are
('sympy','numpy','scipy.sparse'). Currently only 'sympy' is
implemented.
Returns
=======
result : list
A list that consists of primitive states and their probabilities.
Examples
========
>>> from sympy.physics.quantum.qubit import Qubit, measure_all
>>> from sympy.physics.quantum.gate import H, X, Y, Z
>>> from sympy.physics.quantum.qapply import qapply
>>> c = H(0)*H(1)*Qubit('00')
>>> c
H(0)*H(1)*|00>
>>> q = qapply(c)
>>> measure_all(q)
[(|00>, 1/4), (|01>, 1/4), (|10>, 1/4), (|11>, 1/4)]
"""
m = qubit_to_matrix(qubit, format)
if format == 'sympy':
results = []
if normalize:
m = m.normalized()
size = max(m.shape) # Max of shape to account for bra or ket
nqubits = int(math.log(size)/math.log(2))
for i in range(size):
if m[i] != 0.0:
results.append(
(Qubit(IntQubit(i, nqubits=nqubits)), m[i]*conjugate(m[i]))
)
return results
else:
raise NotImplementedError(
"This function can't handle non-sympy matrix formats yet"
)
[docs]def measure_partial(qubit, bits, format='sympy', normalize=True):
"""Perform a partial ensemble measure on the specified qubits.
Parameters
==========
qubits : Qubit
The qubit to measure. This can be any Qubit or a linear combination
of them.
bits : tuple
The qubits to measure.
format : str
The format of the intermediate matrices to use. Possible values are
('sympy','numpy','scipy.sparse'). Currently only 'sympy' is
implemented.
Returns
=======
result : list
A list that consists of primitive states and their probabilities.
Examples
========
>>> from sympy.physics.quantum.qubit import Qubit, measure_partial
>>> from sympy.physics.quantum.gate import H, X, Y, Z
>>> from sympy.physics.quantum.qapply import qapply
>>> c = H(0)*H(1)*Qubit('00')
>>> c
H(0)*H(1)*|00>
>>> q = qapply(c)
>>> measure_partial(q, (0,))
[(sqrt(2)*|00>/2 + sqrt(2)*|10>/2, 1/2), (sqrt(2)*|01>/2 + sqrt(2)*|11>/2, 1/2)]
"""
m = qubit_to_matrix(qubit, format)
if isinstance(bits, (SYMPY_INTS, Integer)):
bits = (int(bits),)
if format == 'sympy':
if normalize:
m = m.normalized()
possible_outcomes = _get_possible_outcomes(m, bits)
# Form output from function.
output = []
for outcome in possible_outcomes:
# Calculate probability of finding the specified bits with
# given values.
prob_of_outcome = 0
prob_of_outcome += (outcome.H*outcome)[0]
# If the output has a chance, append it to output with found
# probability.
if prob_of_outcome != 0:
if normalize:
next_matrix = matrix_to_qubit(outcome.normalized())
else:
next_matrix = matrix_to_qubit(outcome)
output.append((
next_matrix,
prob_of_outcome
))
return output
else:
raise NotImplementedError(
"This function can't handle non-sympy matrix formats yet"
)
[docs]def measure_partial_oneshot(qubit, bits, format='sympy'):
"""Perform a partial oneshot measurement on the specified qubits.
A oneshot measurement is equivalent to performing a measurement on a
quantum system. This type of measurement does not return the probabilities
like an ensemble measurement does, but rather returns *one* of the
possible resulting states. The exact state that is returned is determined
by picking a state randomly according to the ensemble probabilities.
Parameters
----------
qubits : Qubit
The qubit to measure. This can be any Qubit or a linear combination
of them.
bits : tuple
The qubits to measure.
format : str
The format of the intermediate matrices to use. Possible values are
('sympy','numpy','scipy.sparse'). Currently only 'sympy' is
implemented.
Returns
-------
result : Qubit
The qubit that the system collapsed to upon measurement.
"""
import random
m = qubit_to_matrix(qubit, format)
if format == 'sympy':
m = m.normalized()
possible_outcomes = _get_possible_outcomes(m, bits)
# Form output from function
random_number = random.random()
total_prob = 0
for outcome in possible_outcomes:
# Calculate probability of finding the specified bits
# with given values
total_prob += (outcome.H*outcome)[0]
if total_prob >= random_number:
return matrix_to_qubit(outcome.normalized())
else:
raise NotImplementedError(
"This function can't handle non-sympy matrix formats yet"
)
def _get_possible_outcomes(m, bits):
"""Get the possible states that can be produced in a measurement.
Parameters
----------
m : Matrix
The matrix representing the state of the system.
bits : tuple, list
Which bits will be measured.
Returns
-------
result : list
The list of possible states which can occur given this measurement.
These are un-normalized so we can derive the probability of finding
this state by taking the inner product with itself
"""
# This is filled with loads of dirty binary tricks...You have been warned
size = max(m.shape) # Max of shape to account for bra or ket
nqubits = int(math.log(size, 2) + .1) # Number of qubits possible
# Make the output states and put in output_matrices, nothing in them now.
# Each state will represent a possible outcome of the measurement
# Thus, output_matrices[0] is the matrix which we get when all measured
# bits return 0. and output_matrices[1] is the matrix for only the 0th
# bit being true
output_matrices = []
for i in range(1 << len(bits)):
output_matrices.append(zeros(2**nqubits, 1))
# Bitmasks will help sort how to determine possible outcomes.
# When the bit mask is and-ed with a matrix-index,
# it will determine which state that index belongs to
bit_masks = []
for bit in bits:
bit_masks.append(1 << bit)
# Make possible outcome states
for i in range(2**nqubits):
trueness = 0 # This tells us to which output_matrix this value belongs
# Find trueness
for j in range(len(bit_masks)):
if i & bit_masks[j]:
trueness += j + 1
# Put the value in the correct output matrix
output_matrices[trueness][i] = m[i]
return output_matrices
[docs]def measure_all_oneshot(qubit, format='sympy'):
"""Perform a oneshot ensemble measurement on all qubits.
A oneshot measurement is equivalent to performing a measurement on a
quantum system. This type of measurement does not return the probabilities
like an ensemble measurement does, but rather returns *one* of the
possible resulting states. The exact state that is returned is determined
by picking a state randomly according to the ensemble probabilities.
Parameters
----------
qubits : Qubit
The qubit to measure. This can be any Qubit or a linear combination
of them.
format : str
The format of the intermediate matrices to use. Possible values are
('sympy','numpy','scipy.sparse'). Currently only 'sympy' is
implemented.
Returns
-------
result : Qubit
The qubit that the system collapsed to upon measurement.
"""
import random
m = qubit_to_matrix(qubit)
if format == 'sympy':
m = m.normalized()
random_number = random.random()
total = 0
result = 0
for i in m:
total += i*i.conjugate()
if total > random_number:
break
result += 1
return Qubit(IntQubit(result, int(math.log(max(m.shape), 2) + .1)))
else:
raise NotImplementedError(
"This function can't handle non-sympy matrix formats yet"
)
