1568 lines
53 KiB
Python
1568 lines
53 KiB
Python
import random
|
||
from collections import defaultdict
|
||
from functools import reduce
|
||
from numbers import Complex
|
||
from typing import Union
|
||
|
||
import numpy as np
|
||
import sympy as sp
|
||
|
||
from load_test import sizeof_fmt
|
||
|
||
ListOrNdarray = Union[list, np.ndarray]
|
||
|
||
REPR_EMPTY_SET = "Ø"
|
||
REPR_TENSOR_OP = "⊗"
|
||
REPR_TARGET = "⊕"
|
||
REPR_GREEK_BETA = "β"
|
||
REPR_GREEK_PSI = "ψ"
|
||
REPR_GREEK_PHI = "φ"
|
||
REPR_MATH_KET = "⟩"
|
||
REPR_MATH_BRA = "⟨"
|
||
REPR_MATH_SQRT = "√"
|
||
REPR_MATH_SUBSCRIPT_NUMBERS = "₀₁₂₃₄₅₆₇₈₉"
|
||
|
||
# Keep a reference to already used states for naming
|
||
UNIVERSE_STATES = []
|
||
|
||
|
||
class Matrix(object):
|
||
"""Wraps a Matrix... it's for my understanding, this could easily
|
||
probably be np.array"""
|
||
|
||
def __init__(self, m: ListOrNdarray = None, *args, **kwargs):
|
||
"""
|
||
Can be initialized with a matrix, e.g. for |0> this is [[0],[1]]
|
||
:param m: a matrix representing the quantum state
|
||
:param name: the name of the state (optional)
|
||
"""
|
||
if m is None:
|
||
self.m = np.array([])
|
||
elif type(m) is list:
|
||
self.m = np.array(m)
|
||
elif type(m) is np.ndarray:
|
||
self.m = m
|
||
else:
|
||
raise TypeError("m needs to be a list or ndarray type")
|
||
|
||
def __add__(self, other):
|
||
return Matrix(self.m + other.m)
|
||
|
||
def __sub__(self, other):
|
||
return Matrix(self.m - other.m)
|
||
|
||
def __eq__(self, other):
|
||
if isinstance(other, Complex):
|
||
return bool(np.allclose(self.m, other))
|
||
return np.allclose(self.m, other.m)
|
||
|
||
def __mul_or_kron__(self, other):
|
||
"""
|
||
Multiplication with a number is Linear op;
|
||
with another state is a composition via the Kronecker/Tensor product
|
||
"""
|
||
if isinstance(other, Complex):
|
||
return Matrix(self.m * other)
|
||
elif isinstance(other, Matrix):
|
||
return Matrix(np.kron(self.m, other.m))
|
||
raise NotImplementedError
|
||
|
||
def __rmul__(self, other):
|
||
return self.__mul_or_kron__(other)
|
||
|
||
def __mul__(self, other):
|
||
return self.__mul_or_kron__(other)
|
||
|
||
def __len__(self):
|
||
return len(self.m)
|
||
|
||
def inner(self, other):
|
||
"""Define inner product - a.k.a dot product, scalar product - <0|0>"""
|
||
return Matrix(np.dot(self._conjugate_transpose(), other.m))
|
||
|
||
def dot(self, other):
|
||
"""Alias to inner"""
|
||
return self.inner(other)
|
||
|
||
def scalar(self, other):
|
||
"""Alias to inner"""
|
||
return self.inner(other)
|
||
|
||
def outer(self, other):
|
||
"""Define outer product |0><0|
|
||
https://en.wikipedia.org/wiki/Outer_product
|
||
|
||
Given two vectors u and v, their outer product u ⊗ v is defined
|
||
as the m × n matrix A obtained by multiplying each element of u
|
||
by each element of v
|
||
|
||
The outer product u ⊗ v is equivalent to a matrix multiplication uvT,
|
||
provided that u is represented as a m × 1 column vector and v as a
|
||
n × 1 column vector (which makes vT a row vector).
|
||
"""
|
||
return Matrix(np.outer(self.m, other.m))
|
||
|
||
def x(self, other):
|
||
"""Define outer product |0><0| looks like |0x0| which is 0.x(0)"""
|
||
return self.outer(other)
|
||
|
||
def kron(self, other):
|
||
"""Define kronicker product - |0> ⊗ |0> = |00>"""
|
||
return Matrix(np.kron(self.m, other.m))
|
||
|
||
def _conjugate_transpose(self):
|
||
return self.m.transpose().conjugate()
|
||
|
||
def _complex_conjugate(self):
|
||
return self.m.conjugate()
|
||
|
||
def conjugate_transpose(self):
|
||
return Matrix(self._conjugate_transpose())
|
||
|
||
def complex_conjugate(self):
|
||
return Matrix(self._complex_conjugate())
|
||
|
||
@property
|
||
def human_m(self):
|
||
return humanize(self.m)
|
||
|
||
def __repr__(self):
|
||
return str(self.m)
|
||
|
||
|
||
MatrixOrList = Union[Matrix, ListOrNdarray]
|
||
|
||
|
||
class HorizontalVector(Matrix):
|
||
"""Horizontal vector is basically a list"""
|
||
|
||
def __init__(self, m: MatrixOrList = None, *args, **kwargs):
|
||
super().__init__(m, *args, **kwargs)
|
||
if type(m) in [Matrix]:
|
||
super().__init__(m.m, *args, **kwargs)
|
||
else:
|
||
super().__init__(m, *args, **kwargs)
|
||
if not self._is_horizontal_vector():
|
||
raise TypeError("Not a horizontal vector")
|
||
|
||
def _is_horizontal_vector(self):
|
||
return len(self.m.shape) == 1
|
||
|
||
|
||
class Vector(Matrix):
|
||
def __init__(self, m: MatrixOrList = None, *args, **kwargs):
|
||
super().__init__(m, *args, **kwargs)
|
||
if type(m) in [Matrix]:
|
||
super().__init__(m.m, *args, **kwargs)
|
||
else:
|
||
super().__init__(m, *args, **kwargs)
|
||
if not self._is_vector():
|
||
raise TypeError("Not a vector")
|
||
|
||
def _is_vector(self):
|
||
return self.m.shape[1] == 1
|
||
|
||
|
||
VectorOrList = Union[Vector, MatrixOrList]
|
||
|
||
|
||
class State(Vector):
|
||
def __init__(self, m: VectorOrList = None, name: str = '',
|
||
allow_unnormalized=False, *args, **kwargs):
|
||
"""State vector representing quantum state"""
|
||
if type(m) in [Vector, Matrix, State]:
|
||
super().__init__(m.m, *args, **kwargs)
|
||
else:
|
||
super().__init__(m, *args, **kwargs)
|
||
self.name = name
|
||
self.measurement_result = None
|
||
self.last_basis = None
|
||
if not allow_unnormalized and not self._is_normalized():
|
||
raise TypeError("Not a normalized state vector")
|
||
|
||
def __hash__(self):
|
||
return hash(str(self.m))
|
||
|
||
def __le__(self, other):
|
||
return repr(self) < repr(other)
|
||
|
||
@staticmethod
|
||
def normalize(vector: Vector):
|
||
"""Normalize a state by dividing by the square root of sum of the
|
||
squares of states"""
|
||
norm_coef = np.sqrt(np.sum(np.abs(vector.m ** 2)))
|
||
if norm_coef == 0:
|
||
raise TypeError("zero-sum vector")
|
||
new_m = vector.m / norm_coef
|
||
return State(new_m)
|
||
|
||
def _is_normalized(self):
|
||
return np.isclose(np.sum(np.abs(self.m ** 2)), 1.0)
|
||
|
||
@classmethod
|
||
def from_bloch_angles(cls, theta, phi):
|
||
"""Creates a state from angles in Bloch sphere"""
|
||
m0 = np.cos(theta / 2)
|
||
m1 = np.sin(theta / 2) * np.power(np.e, (1j * phi))
|
||
m = m0 * _0 + m1 * _1
|
||
return cls(m.m)
|
||
|
||
def to_bloch_angles(self):
|
||
"""Returns the angles of this state on the Bloch sphere"""
|
||
if not self.m.shape == (2, 1):
|
||
raise Exception("State needs to describe only 1 qubit on the bloch sphere (2x1 matrix)")
|
||
m0, m1 = self.m[0][0], self.m[1][0]
|
||
|
||
# theta is between 0 and pi
|
||
theta = 2 * np.arccos(m0)
|
||
if theta < 0:
|
||
theta += np.pi
|
||
assert 0 <= theta <= np.pi
|
||
|
||
div = np.sin(theta / 2)
|
||
if div == 0:
|
||
# here is doesn't matter what phi is as phase at the poles is arbitrary
|
||
phi = 0
|
||
else:
|
||
exp = m1 / div
|
||
phi = np.log(complex(exp)) / 1j
|
||
if phi < 0:
|
||
phi += 2 * np.pi
|
||
assert 0 <= phi <= 2 * np.pi
|
||
return theta, phi
|
||
|
||
def to_ket(self):
|
||
return self.conjugate_transpose()
|
||
|
||
def rotate_x(self, theta):
|
||
return Rx(theta).on(self)
|
||
|
||
def rotate_y(self, theta):
|
||
return Ry(theta).on(self)
|
||
|
||
def rotate_z(self, theta):
|
||
return Rz(theta).on(self)
|
||
|
||
def __repr__(self):
|
||
if not self.name:
|
||
if np.count_nonzero(self.m == 1):
|
||
fmt_str = self.get_fmt_of_element()
|
||
index = np.where(self.m == [1])[0][0]
|
||
self.name = fmt_str.format(index)
|
||
else:
|
||
for well_known_state in well_known_states:
|
||
try:
|
||
if self == well_known_state:
|
||
return repr(well_known_state)
|
||
except:
|
||
...
|
||
for used_state in UNIVERSE_STATES:
|
||
try:
|
||
if self == used_state:
|
||
return repr(used_state)
|
||
except:
|
||
...
|
||
next_state_sub = ''.join(
|
||
[REPR_MATH_SUBSCRIPT_NUMBERS[int(d)] for d in
|
||
str(len(UNIVERSE_STATES))])
|
||
self.name = '{}{}'.format(REPR_GREEK_PSI, next_state_sub)
|
||
UNIVERSE_STATES.append(self)
|
||
state_name = '|{}>'.format(self.name)
|
||
return state_name
|
||
|
||
def norm(self):
|
||
"""Norm/Length of the vector = sqrt(<self|self>)"""
|
||
return self.length()
|
||
|
||
def length(self):
|
||
"""Norm/Length of the vector = sqrt(<self|self>)"""
|
||
return np.sqrt((self.inner(self)).m).item(0)
|
||
|
||
def is_orthogonal(self, other):
|
||
"""If the inner (dot) product is zero, this vector is orthogonal to
|
||
other"""
|
||
return self.inner(other) == 0
|
||
|
||
def get_prob_from_measurement_op(self, m_op):
|
||
assert type(m_op) is MeasurementOperator
|
||
return m_op.get_prob(self)
|
||
|
||
def get_computational_basis_vector(self, j):
|
||
return Vector([[1] if i == int(j) else [0] for i in range(len(self.m))])
|
||
|
||
def get_prob(self, j=None, basis=None):
|
||
"""pr(j) = |<e_j|self>|^2
|
||
"""
|
||
if basis is None:
|
||
assert j is not None
|
||
basis = self.get_computational_basis_vector(j)
|
||
m = MeasurementOperator.create_from_basis(basis)
|
||
return m.get_prob(self)
|
||
|
||
def get_fmt_of_element(self):
|
||
return "{:0" + str(int(np.ceil(np.log2(len(self))))) + "b}"
|
||
|
||
def measure(self, basis=None, allow_empty=False):
|
||
"""
|
||
m_ops: a set of MeasurementOperators.
|
||
e.g. for computational basis, the m_ops is [|0><0|, |1><1|]
|
||
Measures in the computational basis.
|
||
Irreversable operation. Measuring again will result in the same result
|
||
TODO: Should we memoize per basis?
|
||
Currently it's memoized if the basis is the same
|
||
If it's measured twice, should it return the same state?
|
||
What about if measured twice but in different bases?
|
||
E.g. measure1 -> computation -> A
|
||
measure2 -> +/- basis -> B
|
||
measure3 -> computation -> should it return A again or
|
||
random weighted?
|
||
measure4 -> +/- basis -> should it return B again or
|
||
random weighted?
|
||
:return: binary representation of the measured qubit (e.g. "011")
|
||
"""
|
||
weights, m_ops = [], []
|
||
if not basis:
|
||
basis = []
|
||
for j in range(len(self)):
|
||
# Default is computational basis
|
||
e_j = self.get_computational_basis_vector(j)
|
||
m = MeasurementOperator.create_from_basis(e_j)
|
||
m_ops.append(m)
|
||
basis.append(State(e_j))
|
||
else:
|
||
assert len(basis) == len(self)
|
||
for i, e_j in enumerate(basis):
|
||
# All should be orthogonal
|
||
if i != 0:
|
||
assert e_j.is_orthogonal(basis[0])
|
||
m_ops.append(MeasurementOperator.create_from_basis(e_j))
|
||
|
||
if self.measurement_result and self.last_basis == basis:
|
||
return self.measurement_result
|
||
|
||
for m_op in m_ops:
|
||
prob = self.get_prob_from_measurement_op(m_op)
|
||
weights += [prob]
|
||
|
||
empty_choices, empty_weights = [], []
|
||
if allow_empty:
|
||
empty_prob = 1.0 - sum(weights)
|
||
empty_weights = [empty_prob]
|
||
empty_choices = [REPR_EMPTY_SET]
|
||
else:
|
||
# TODO: This may be wrong
|
||
# normalize the weights
|
||
weights_sum = np.sum(weights)
|
||
if not np.isclose(weights_sum, 1.0):
|
||
weights = list(np.array(weights) / weights_sum)
|
||
assert np.isclose(np.sum(weights), 1.0)
|
||
|
||
weights = empty_weights + weights
|
||
self.measurement_result = random.choices(empty_choices + basis, weights)[0]
|
||
self.last_basis = basis
|
||
return self.measurement_result
|
||
|
||
def measure_with_op(self, mo):
|
||
"""
|
||
Measures with a measurement operator mo
|
||
TODO: Can't define `mo: MeasurementOperator` because in python you
|
||
can't declare classes before defining them
|
||
"""
|
||
m = mo.on(self).m / np.sqrt(mo.get_prob(self))
|
||
return State(m)
|
||
|
||
def measure_partial(self, qubit_n):
|
||
"""Partial measurement of state
|
||
Measures the n-th qubit with probability sum(a_n),
|
||
adjusting the probability of the rest of the state
|
||
"""
|
||
max_qubits = int(np.log2(len(self)))
|
||
if not (0 < qubit_n <= max_qubits):
|
||
raise Exception("Partial measurement of qubit_n must be between 1 and {}".format(max_qubits))
|
||
format_str = self.get_fmt_of_element()
|
||
# e.g. for state |000>:
|
||
# ['000', '001', '010', '011', '100', '101', '110', '111']
|
||
bin_repr = [format_str.format(i) for i in range(len(self))]
|
||
# e.g. for qbit_n = 2
|
||
# [0, 0, 1, 1, 0, 0, 1, 1]
|
||
partial_measurement_of_qbit = [int(b[qubit_n - 1]) for b in bin_repr]
|
||
# [0, 1, 4, 5]
|
||
weights, choices = defaultdict(list), defaultdict(list)
|
||
for result in [1, 0]:
|
||
indexes_for_p_0 = [i for i, index in enumerate(partial_measurement_of_qbit) if index == result]
|
||
weights[result] = [self.get_prob(j) for j in indexes_for_p_0]
|
||
choices[result] = [format_str.format(i) for i in indexes_for_p_0]
|
||
weights_01 = [sum(weights[0]), sum(weights[1])]
|
||
measurement_result = random.choices([0, 1], weights_01)[0]
|
||
normalization_factor = np.sqrt(sum([np.abs(i) ** 2 for i in weights[measurement_result]]))
|
||
new_m = weights[measurement_result] / normalization_factor
|
||
return State(s('|' + str(measurement_result) + '>')), State(new_m.reshape((len(new_m), 1)))
|
||
|
||
def pretty_print(self):
|
||
format_str = self.get_fmt_of_element() + " | {}"
|
||
|
||
for i, element in enumerate(self.m):
|
||
print(format_str.format(i, humanize_num(element[0])))
|
||
|
||
@classmethod
|
||
def from_string(cls, q):
|
||
try:
|
||
bin_repr = q[1:-1]
|
||
dec_bin = int(bin_repr, 2)
|
||
bin_length = 2 ** len(bin_repr)
|
||
arr = [[1] if i == dec_bin else [0] for i in range(bin_length)]
|
||
except:
|
||
raise Exception(
|
||
"State from string should be of the form |00..01> with all "
|
||
"numbers either 0 or 1")
|
||
return cls(arr)
|
||
|
||
@staticmethod
|
||
def get_random_state_angles():
|
||
phi = np.random.uniform(0, 2 * np.pi)
|
||
t = np.random.uniform(0, 1)
|
||
theta = np.arccos(1 - 2 * t)
|
||
return theta, phi
|
||
|
||
def get_bloch_coordinates(self):
|
||
theta, phi = self.to_bloch_angles()
|
||
x = np.sin(theta) * np.cos(phi)
|
||
y = np.sin(theta) * np.sin(phi)
|
||
z = np.cos(theta)
|
||
return [x, y, z]
|
||
|
||
|
||
class Ket(State):
|
||
def __init__(self, *args, **kwargs):
|
||
"""Alias: Ket is a state"""
|
||
State.__init__(self, *args, **kwargs)
|
||
|
||
|
||
def s(q, name=None):
|
||
"""Helper method for creating state easily"""
|
||
if type(q) == str:
|
||
if q[0] == '|' and q[-1] == '>':
|
||
# e.g. |000>
|
||
state = State.from_string(q)
|
||
state.name = name
|
||
return state
|
||
elif q[0] == '|' and q[-1] == '|' and "><" in q:
|
||
# e.g. |0><0|
|
||
f_qbit, s_qbit = q.split("<")[0], q.split(">")[1]
|
||
if len(f_qbit) == len(s_qbit):
|
||
s_qbit = '|' + s_qbit[1:-1] + '>'
|
||
fs, ss = State.from_string(f_qbit), State.from_string(s_qbit)
|
||
return DensityMatrix(fs.x(ss).m, name=name)
|
||
elif type(q) == list:
|
||
# e.g. s([1,0]) => |0>
|
||
return State(np.reshape(q, (len(q), 1)), name=name)
|
||
return State(q, name=name)
|
||
|
||
|
||
def humanize_num(fl, tolerance=1e-3):
|
||
if np.abs(fl.imag) < tolerance:
|
||
fl = fl.real
|
||
if np.abs(fl.real) < tolerance:
|
||
fl = 0 + 1j * fl.imag
|
||
try:
|
||
return sp.nsimplify(fl, [sp.pi], tolerance, full=True)
|
||
except:
|
||
return sp.nsimplify(fl, [sp.pi], tolerance)
|
||
|
||
|
||
def pp(m, tolerance=1e-3):
|
||
for element in m:
|
||
print(humanize_num(element, tolerance=tolerance))
|
||
|
||
|
||
def humanize(m):
|
||
rv = []
|
||
for element in m:
|
||
if type(element) in [np.ndarray, list]:
|
||
rv.append(humanize(element))
|
||
else:
|
||
rv.append(humanize_num(element))
|
||
return rv
|
||
|
||
|
||
class Operator(object):
|
||
def __init__(self, func: sp.Lambda, *args, **kwargs):
|
||
"""An Operator turns one function into another"""
|
||
self.func = func
|
||
|
||
def on(self, *args, **kwargs):
|
||
return self(*args, **kwargs)
|
||
|
||
def __call__(self, *args, **kwargs):
|
||
return self.func(*args, **kwargs)
|
||
|
||
|
||
class LinearOperator(Operator):
|
||
def __init__(self, func: sp.Lambda, *args, **kwargs):
|
||
"""Linear operators satisfy f(x+y) = f(x) + f(y) and a*f(x) = f(a*x)"""
|
||
super().__init__(func, *args, **kwargs)
|
||
if not self._is_linear():
|
||
raise TypeError("Not a linear operator")
|
||
|
||
def _is_linear(self):
|
||
# A function is (jointly) linear in a given set of variables
|
||
# if all second-order derivatives are identically zero (including
|
||
# mixed ones).
|
||
# https://stackoverflow.com/questions/36283548/check-if-an-equation
|
||
# -is-linear-for-a-specific-set-of-variables
|
||
expr, vars_ = self.func.expr, self.func.variables
|
||
for x in vars_:
|
||
for y in vars_:
|
||
try:
|
||
if not sp.Eq(sp.diff(expr, x, y), 0):
|
||
return False
|
||
except TypeError:
|
||
return False
|
||
return True
|
||
|
||
|
||
def test_linear_operator():
|
||
# Operators turn one vector into another
|
||
# the times 2 operator should return the times two multiplication
|
||
_x = sp.Symbol('x')
|
||
_times_2 = LinearOperator(sp.Lambda(_x, 2 * _x))
|
||
assert _times_2.on(5) == 10
|
||
assert _times_2(5) == 10
|
||
|
||
assert_raises(TypeError, "Not a linear operator", LinearOperator,
|
||
sp.Lambda(_x, _x ** 2))
|
||
|
||
|
||
class SquareMatrix(Matrix):
|
||
def __init__(self, m: ListOrNdarray, *args, **kwargs):
|
||
super().__init__(m, *args, **kwargs)
|
||
if not self._is_square():
|
||
raise TypeError("Not a Square matrix")
|
||
|
||
def _is_square(self):
|
||
return self.m.shape[0] == self.m.shape[1]
|
||
|
||
|
||
class LinearTransformation(LinearOperator, Matrix):
|
||
def __init__(self, m: ListOrNdarray, func=None, name: str = '', *args,
|
||
**kwargs):
|
||
"""LinearTransformation (or linear map) inherits from both
|
||
LinearOperator and a Matrix
|
||
It is used to act on a State vector by defining the operator to be
|
||
the dot product"""
|
||
self.name = name
|
||
Matrix.__init__(self, m=m, *args, **kwargs)
|
||
self.func = func or self.operator_func
|
||
LinearOperator.__init__(self, func=func, *args, **kwargs)
|
||
|
||
def _is_linear(self):
|
||
# Every matrix transformation is a linear transformation
|
||
# https://www.mathbootcamps.com/proof-every-matrix-transformation-is
|
||
# -a-linear-transformation/
|
||
return True
|
||
|
||
def operator_func(self, other):
|
||
return Vector(np.dot(self.m, other.m))
|
||
|
||
def get_eigens(self):
|
||
""" Returns (eigenvalue, eigenvector)
|
||
M|v> = λ|v> ->
|
||
M :is the operator (self)
|
||
|v> :is eigenstate of M
|
||
λ :is the corresponding eigenvalue
|
||
"""
|
||
eigenvalues, eigenvectors = np.linalg.eig(self.m)
|
||
rv = []
|
||
for i in range(0, len(eigenvectors)):
|
||
eigenvalue = eigenvalues[i]
|
||
eigenvector = HorizontalVector(eigenvectors[:, i])
|
||
rv.append((eigenvalue, eigenvector))
|
||
return rv
|
||
|
||
|
||
class UnitaryMatrix(SquareMatrix):
|
||
def __init__(self, m: ListOrNdarray, *args, **kwargs):
|
||
"""Represents a Unitary matrix that satisfies UU+ = I"""
|
||
super().__init__(m, *args, **kwargs)
|
||
if not self._is_unitary():
|
||
raise TypeError("Not a Unitary matrix")
|
||
|
||
def _is_unitary(self):
|
||
"""Checks if the matrix product of itself with conjugate transpose is
|
||
Identity UU+ = I"""
|
||
UU_ = np.dot(self._conjugate_transpose(), self.m)
|
||
I = np.eye(self.m.shape[0])
|
||
return np.isclose(UU_, I).all()
|
||
|
||
|
||
class HermitianMatrix(SquareMatrix):
|
||
def __init__(self, m: ListOrNdarray, *args, **kwargs):
|
||
"""Represents a Hermitian matrix that satisfies U=U+"""
|
||
super().__init__(m, *args, **kwargs)
|
||
if not self._is_hermitian():
|
||
raise TypeError("Not a Hermitian matrix")
|
||
|
||
def _is_hermitian(self):
|
||
"""Checks if its equal to the conjugate transpose U = U+"""
|
||
return np.isclose(self.m, self._conjugate_transpose()).all()
|
||
|
||
|
||
class UnitaryOperator(LinearTransformation, UnitaryMatrix):
|
||
def __init__(self, m: ListOrNdarray, name: str = '', partials=None,
|
||
decomposition=None, *args, **kwargs):
|
||
"""UnitaryOperator inherits from both LinearTransformation and a
|
||
Unitary matrix
|
||
It is used to act on a State vector by defining the operator to be
|
||
the dot product
|
||
|
||
Partials - a list of partial operators that are used in a multi-qubit
|
||
UnitaryOperator to define the operator on each qubit. Each element of
|
||
the list is the Nth partial that is used, i.e. for the first - |0><0|,
|
||
for the second - |1><1|
|
||
"""
|
||
if np.shape(m) != (2, 2) and partials is None and decomposition is None:
|
||
raise Exception("Please define partials or decomposition in "
|
||
"a non-single operator")
|
||
UnitaryMatrix.__init__(self, m=m, *args, **kwargs)
|
||
LinearTransformation.__init__(self, m=m, func=self.operator_func, *args,
|
||
**kwargs)
|
||
self.name = name
|
||
self.partials = partials
|
||
self.decomposition = decomposition
|
||
|
||
def operator_func(self, other, which_qbit=None):
|
||
this_cols, other_rows = np.shape(self.m)[1], np.shape(other.m)[0]
|
||
if this_cols == other_rows and which_qbit is None:
|
||
return State(np.dot(self.m, other.m))
|
||
if which_qbit is None:
|
||
raise Exception("Operating dim-{} operator on a dim-{} state. "
|
||
"Please specify which_qubit to operate on".format(
|
||
this_cols, other_rows))
|
||
total_qbits = int(np.log2(other_rows))
|
||
if type(which_qbit) is list and len(which_qbit) == 1:
|
||
which_qbit = which_qbit[0]
|
||
if type(which_qbit) is int:
|
||
# single qubit-gate
|
||
assert this_cols == 2
|
||
assert 0 <= which_qbit < total_qbits
|
||
extended_m = [self if i == which_qbit else I for i in
|
||
range(total_qbits)]
|
||
new_m = np.prod(extended_m).m
|
||
elif type(which_qbit) is list:
|
||
if self.decomposition:
|
||
return self.decomposition(other, *which_qbit)
|
||
else:
|
||
# single or multiple qubit-gate
|
||
assert 1 <= len(which_qbit) <= total_qbits
|
||
assert len(which_qbit) == len(self.partials)
|
||
assert all([q < total_qbits for q in which_qbit])
|
||
this_gate_len = 2 ** (len(self.partials) - 1)
|
||
bins = generate_bins(this_gate_len)
|
||
extended_m, next_control = [[] for _ in bins], 0
|
||
partial_mapping = dict(zip(which_qbit, self.partials))
|
||
for qbit in range(total_qbits):
|
||
if qbit not in partial_mapping:
|
||
for i in range(this_gate_len):
|
||
extended_m[i].append(I)
|
||
else:
|
||
this_partial = partial_mapping[qbit]
|
||
# TODO: This works only with C_partial :(((
|
||
if this_partial == C_partial:
|
||
for i in range(this_gate_len):
|
||
bin_dig = bins[i][next_control]
|
||
extended_m[i].append(
|
||
s("|" + bin_dig + "><" + bin_dig + "|"))
|
||
next_control += 1
|
||
else:
|
||
for i in range(this_gate_len - 1):
|
||
extended_m[i].append(I)
|
||
extended_m[-1].append(this_partial)
|
||
new_m = sum([np.prod(e).m for e in extended_m])
|
||
else:
|
||
raise Exception(
|
||
"which_qubit needs to be either an int of N-th qubit or list")
|
||
|
||
return State(np.dot(new_m, other.m))
|
||
|
||
def __repr__(self):
|
||
if self.name:
|
||
return '-{}-'.format(self.name)
|
||
return str(self.m)
|
||
|
||
|
||
class Gate(UnitaryOperator):
|
||
def __init__(self, *args, **kwargs):
|
||
"""Alias: Gates are Unitary Operators"""
|
||
UnitaryOperator.__init__(self, *args, **kwargs)
|
||
|
||
|
||
class HermitianOperator(LinearTransformation, HermitianMatrix):
|
||
def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs):
|
||
"""HermitianOperator inherits from both LinearTransformation and a
|
||
Hermitian matrix
|
||
It is used to act on a State vector by defining the operator to be
|
||
the dot product"""
|
||
self.name = name
|
||
HermitianMatrix.__init__(self, m=m, *args, **kwargs)
|
||
LinearOperator.__init__(self, func=self.operator_func, *args, **kwargs)
|
||
|
||
def operator_func(self, other):
|
||
"""This might not return a normalized state vector so don't wrap it
|
||
in State"""
|
||
return Vector(np.dot(self.m, other.m))
|
||
|
||
def __repr__(self):
|
||
if self.name:
|
||
return '-{}-'.format(self.name)
|
||
return str(self.m)
|
||
|
||
|
||
class DensityMatrix(HermitianOperator):
|
||
def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs):
|
||
"""The density matrix is a representation of a linear operator called
|
||
the density operator. The density matrix is obtained from the density
|
||
operator by choice of basis in the underlying space. In practice,
|
||
the terms density matrix and density operator are often used
|
||
interchangeably. Both matrix and operator are self-adjoint
|
||
(or Hermitian), positive semi-definite, of trace one, and may
|
||
be infinite-dimensional."""
|
||
self.name = name
|
||
HermitianOperator.__init__(self, m=m, *args, **kwargs)
|
||
if not self._is_density():
|
||
raise TypeError("Not a Density matrix")
|
||
|
||
def _is_density(self):
|
||
"""Checks if it trace is 1"""
|
||
return np.isclose(np.trace(self.m), 1.0)
|
||
|
||
|
||
"""
|
||
Define States and Operators
|
||
"""
|
||
|
||
# EMPTY STATE
|
||
_E = Vector([[0],
|
||
[0]],
|
||
name=REPR_EMPTY_SET)
|
||
|
||
_0 = State([[1],
|
||
[0]],
|
||
name='0')
|
||
|
||
_1 = State([[0],
|
||
[1]],
|
||
name='1')
|
||
|
||
# |+> and |-> states
|
||
_p = State([[1 / np.sqrt(2)],
|
||
[1 / np.sqrt(2)]],
|
||
name='+')
|
||
|
||
_m = State([[1 / np.sqrt(2)],
|
||
[- 1 / np.sqrt(2)]],
|
||
name='-')
|
||
|
||
# 4 Bell states
|
||
b_phi_p = State([[1 / np.sqrt(2)],
|
||
[0],
|
||
[0],
|
||
[1 / np.sqrt(2)]],
|
||
name="{}+".format(REPR_GREEK_PHI))
|
||
|
||
b_psi_p = State([[0],
|
||
[1 / np.sqrt(2)],
|
||
[1 / np.sqrt(2)],
|
||
[0]],
|
||
name="{}+".format(REPR_GREEK_PSI))
|
||
|
||
b_phi_m = State([[1 / np.sqrt(2)],
|
||
[0],
|
||
[0],
|
||
[-1 / np.sqrt(2)]],
|
||
name="{}-".format(REPR_GREEK_PHI))
|
||
|
||
b_psi_m = State([[0],
|
||
[1 / np.sqrt(2)],
|
||
[-1 / np.sqrt(2)],
|
||
[0]],
|
||
name="{}-".format(REPR_GREEK_PSI))
|
||
|
||
# define common bases:
|
||
computational_basis = [_0, _1]
|
||
plus_minus_basis = [_p, _m]
|
||
bell_basis = [b_phi_p, b_psi_p, b_phi_m, b_psi_m]
|
||
|
||
well_known_states = [_p, _m, b_phi_p, b_psi_p, b_phi_m, b_psi_m]
|
||
|
||
_ = I = Gate([[1, 0],
|
||
[0, 1]],
|
||
name="-")
|
||
|
||
X = Gate([[0, 1],
|
||
[1, 0]],
|
||
name="X")
|
||
|
||
Y = Gate([[0, -1j],
|
||
[1j, 0]],
|
||
name="Y")
|
||
|
||
Z = Gate([[1, 0],
|
||
[0, -1]],
|
||
name="Z")
|
||
|
||
# These are rotations that are specified commonly e.g. in
|
||
# https://www.quantum-inspire.com/kbase/rotation-operators/
|
||
# http://www.vcpc.univie.ac.at/~ian/hotlist/qc/talks/bloch-sphere
|
||
# -rotations.pdf
|
||
# and elsewhere DO NOT equate X, Y and Z for theta=np.pi
|
||
# However - they are correct up to a global phase which is all that matters
|
||
# for measurement purposes
|
||
#
|
||
|
||
H = Gate([[1 / np.sqrt(2), 1 / np.sqrt(2)],
|
||
[1 / np.sqrt(2), -1 / np.sqrt(2)], ],
|
||
name="H")
|
||
|
||
# Rotation operators
|
||
# https://www.quantum-inspire.com/kbase/rx-gate/
|
||
Rx = lambda theta: Gate([[np.cos(theta / 2), -1j * np.sin(theta / 2)],
|
||
[-1j * np.sin(theta / 2), np.cos(theta / 2)]],
|
||
name="Rx")
|
||
|
||
Ry = lambda theta: Gate([[np.cos(theta / 2), -np.sin(theta / 2)],
|
||
[np.sin(theta / 2), np.cos(theta / 2)]],
|
||
name="Ry")
|
||
|
||
Rz = lambda theta: Gate([[np.power(np.e, -1j * theta / 2), 0],
|
||
[0, np.power(np.e, 1j * theta / 2)]],
|
||
name="Rz")
|
||
|
||
# See [T-Gate](https://www.quantum-inspire.com/kbase/t-gate/)
|
||
T = lambda phi: Gate([[1, 0],
|
||
[0, np.power(np.e, 1j * phi / 2)]],
|
||
name="T")
|
||
|
||
# How to add a CNOT gate to the Quantum Processor?
|
||
# Imagine if I have to act on a 3-qubit computer and CNOT(q1, q3)
|
||
#
|
||
# Decomposed CNOT :
|
||
# reverse engineered from
|
||
# https://quantumcomputing.stackexchange.com/questions/4252/how-to-derive-the
|
||
# -cnot-matrix-for-a-3-qbit-system-where-the-control-target-qbi
|
||
#
|
||
# CNOT(q1, I, q2):
|
||
# |0><0| x I_2 x I_2 + |1><1| x I_2 x X
|
||
# np.kron(np.kron(np.outer(_0.m, _0.m), np.eye(2)), np.eye(2)) + np.kron(
|
||
# np.kron(np.outer(_1.m, _1.m), np.eye(2)), X.m)
|
||
#
|
||
#
|
||
# CNOT(q1, q2):
|
||
# |0><0| x I + |1><1| x X
|
||
# np.kron(np.outer(_0.m, _0.m), np.eye(2)) + np.kron(np.outer(_1.m, _1.m), X.m)
|
||
# _0.x(_0) * Matrix(I.m) + _1.x(_1) * Matrix(X.m)
|
||
|
||
# Additionally for Toffoli: on a 4qubit system with C-CX (qbits 0,2,3)
|
||
# s("|0><0|") * I * s("|0><0|") * I \
|
||
# + s("|0><0|") * I * s("|1><1|") * I \
|
||
# + s("|1><1|") * I * s("|0><0|") * I \
|
||
# + s("|1><1|") * I * s("|1><1|") * X
|
||
|
||
C_partial = Gate(I.m, name="C")
|
||
x_partial = Gate(X.m, name=REPR_TARGET)
|
||
z_partial = Gate(Z.m, name=REPR_TARGET)
|
||
|
||
CNOT = Gate([
|
||
[1, 0, 0, 0],
|
||
[0, 1, 0, 0],
|
||
[0, 0, 0, 1],
|
||
[0, 0, 1, 0],
|
||
], name="CNOT", partials=[C_partial, x_partial])
|
||
|
||
CZ = Gate([
|
||
[1, 0, 0, 0],
|
||
[0, 1, 0, 0],
|
||
[0, 0, 1, 0],
|
||
[0, 0, 0, -1],
|
||
], name="CZ", partials=[C_partial, z_partial])
|
||
|
||
# Can't do partials for SWAP, but can be decomposed
|
||
# SWAP is decomposed into three CNOTs where the second one is flipped
|
||
SWAP = Gate([
|
||
[1, 0, 0, 0],
|
||
[0, 0, 1, 0],
|
||
[0, 1, 0, 0],
|
||
[0, 0, 0, 1],
|
||
], name="SWAP",
|
||
decomposition=lambda state, x, y:
|
||
CNOT.on(CNOT.on(CNOT.on(state, [x, y]), [y, x]), [x, y])
|
||
)
|
||
|
||
TOFF = Gate([
|
||
[1, 0, 0, 0, 0, 0, 0, 0],
|
||
[0, 1, 0, 0, 0, 0, 0, 0],
|
||
[0, 0, 1, 0, 0, 0, 0, 0],
|
||
[0, 0, 0, 1, 0, 0, 0, 0],
|
||
[0, 0, 0, 0, 1, 0, 0, 0],
|
||
[0, 0, 0, 0, 0, 1, 0, 0],
|
||
[0, 0, 0, 0, 0, 0, 0, 1],
|
||
[0, 0, 0, 0, 0, 0, 1, 0],
|
||
], partials=[C_partial, C_partial, x_partial])
|
||
|
||
|
||
def assert_raises(exception, msg, callable, *args, **kwargs):
|
||
try:
|
||
callable(*args, **kwargs)
|
||
except exception as e:
|
||
assert str(e) == msg
|
||
return
|
||
assert False
|
||
|
||
|
||
def assert_not_raises(exception, msg, callable, *args, **kwargs):
|
||
try:
|
||
callable(*args, **kwargs)
|
||
except exception as e:
|
||
if str(e) == msg:
|
||
assert False
|
||
|
||
|
||
def test_unitary_hermitian():
|
||
# Unitary is UU+ = I; Hermitian is U = U+
|
||
# Matrixes could be either, neither or both
|
||
# Quantum operators (gates) are described *only* by unitary transformations
|
||
# Hermitian operators are used for measurement operators -
|
||
# https://towardsdatascience.com/understanding-basics-of-measurements-in
|
||
# -quantum-computation-4c885879eba0
|
||
h_not_u = [
|
||
[1, 0],
|
||
[0, 2],
|
||
]
|
||
assert_not_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix,
|
||
h_not_u)
|
||
assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, h_not_u)
|
||
|
||
u_not_h = [
|
||
[1, 0],
|
||
[0, 1j],
|
||
]
|
||
assert_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix, u_not_h)
|
||
assert_not_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, u_not_h)
|
||
|
||
u_and_h = [
|
||
[0, 1],
|
||
[1, 0],
|
||
]
|
||
assert_not_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix,
|
||
u_and_h)
|
||
assert_not_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, u_and_h)
|
||
|
||
not_u_not_h = [
|
||
[1, 2],
|
||
[0, 1],
|
||
]
|
||
assert_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix,
|
||
not_u_not_h)
|
||
assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, not_u_not_h)
|
||
|
||
|
||
class MeasurementOperator(HermitianOperator):
|
||
"""Measurement operators are Hermitians: <ψ|M†_m M_m|ψ>"""
|
||
|
||
@classmethod
|
||
def create_from_basis(cls, matrix: Matrix, *args, **kwargs):
|
||
"""returns |M†_m><M_m|"""
|
||
return cls(matrix.conjugate_transpose().x(matrix).m, *args, **kwargs)
|
||
|
||
def get_prob(self, state: State):
|
||
"""Returns result of <ψ|M†_m M_m|ψ>
|
||
"""
|
||
# <ψ|
|
||
state_ct = state.conjugate_transpose()
|
||
# This is: <ψ| . M†_m M_m . |ψ>
|
||
return state_ct.m.dot(self.m.dot(state.m)).item().real
|
||
|
||
|
||
m0 = MeasurementOperator.create_from_basis(_0)
|
||
m1 = MeasurementOperator.create_from_basis(_1)
|
||
|
||
|
||
def test_measurement_ops():
|
||
assert m0 == Matrix([[1, 0],
|
||
[0, 0]])
|
||
assert m1 == Matrix([[0, 0],
|
||
[0, 1]])
|
||
|
||
# p(0) -> probability of measurement to yield a 0
|
||
assert m0.get_prob(_0) == 1.0
|
||
assert m1.get_prob(_0) == 0.0
|
||
assert m0.get_prob(_1) == 0.0
|
||
assert m1.get_prob(_1) == 1.0
|
||
|
||
# Post-state measurement of qubit with operator
|
||
assert _p.measure_with_op(m0) == _0
|
||
assert _p.measure_with_op(m1) == _1
|
||
assert _m.measure_with_op(m0) == _0
|
||
assert _m.measure_with_op(m1) == s([0, -1])
|
||
|
||
|
||
def abs_squared(x):
|
||
return np.abs(x) ** 2
|
||
|
||
|
||
def testRotPauli():
|
||
# From http://www.vcpc.univie.ac.at/~ian/hotlist/qc/talks/bloch-sphere.pdf
|
||
# / slide 11:
|
||
# However, the only measurable quantities are the probabilities
|
||
# |α|**2 and |β|**2, so multiplying the state by an arbitrary factor
|
||
# exp(iγ) (a global phase) has no observable consequences
|
||
assert np.allclose(abs_squared(Rx(np.pi).m), abs_squared(X.m))
|
||
assert np.allclose(abs_squared(Ry(np.pi).m), abs_squared(Y.m))
|
||
assert np.allclose(abs_squared(Rz(np.pi).m), abs_squared(Z.m))
|
||
|
||
|
||
def test_eigenstuff():
|
||
assert LinearTransformation(m=[[1, 0], [0, 0]]).get_eigens() == \
|
||
[(1.0, HorizontalVector([1., 0.])), (0., HorizontalVector([0., 1.]))]
|
||
|
||
|
||
def test_partials():
|
||
# test single qbit state
|
||
assert X.on(s("|000>"), which_qbit=1) == s("|010>")
|
||
assert X.on(s("|000>"), which_qbit=[1]) == s("|010>")
|
||
|
||
# normal 2 qbit state
|
||
assert CNOT.on(s("|00>")) == s("|00>")
|
||
assert CNOT.on(s("|10>")) == s("|11>")
|
||
assert_raises(Exception,
|
||
"Operating dim-4 operator on a dim-8 state. Please specify "
|
||
"which_qubit to operate on",
|
||
CNOT.on, s("|100>"))
|
||
|
||
# Test flipped CNOT
|
||
assert CNOT.on(s("|01>"), which_qbit=[1, 0]) == s("|11>")
|
||
assert CNOT.on(s("|001>"), which_qbit=[2, 0]) == s("|101>")
|
||
|
||
# Test SWAP via 3 successive CNOTs, the second CNOT is flipped
|
||
assert CNOT.on(CNOT.on(CNOT.on(s("|10>")), which_qbit=[1, 0])) == s("|01>")
|
||
|
||
# Test SWAP via 3 successive CNOTs, the 0<->2 are swapped
|
||
assert CNOT.on(CNOT.on(CNOT.on(s("|100>"), [0, 2]), [2, 0]), [0, 2]) == s(
|
||
"|001>")
|
||
|
||
# apply on 0, 1 of 3qbit state
|
||
assert CNOT.on(s("|000>"), which_qbit=[0, 1]) == s("|000>")
|
||
assert CNOT.on(s("|100>"), which_qbit=[0, 1]) == s("|110>")
|
||
|
||
# apply on 1, 2 of 3qbit state
|
||
assert CNOT.on(s("|000>"), which_qbit=[1, 2]) == s("|000>")
|
||
assert CNOT.on(s("|010>"), which_qbit=[1, 2]) == s("|011>")
|
||
|
||
# apply on 0, 2 of 3qbit state
|
||
assert CNOT.on(s("|000>"), which_qbit=[0, 2]) == s("|000>")
|
||
assert CNOT.on(s("|100>"), which_qbit=[0, 2]) == s("|101>")
|
||
|
||
# apply on 0, 2 of 4qbit state
|
||
assert CNOT.on(s("|1000>"), which_qbit=[0, 2]) == s("|1010>")
|
||
|
||
# apply on 0, 3 of 4qbit state
|
||
assert CNOT.on(s("|1000>"), which_qbit=[0, 3]) == s("|1001>")
|
||
|
||
# test SWAP gate
|
||
assert SWAP.on(s("|10>")) == s("|01>")
|
||
assert SWAP.on(s("|01>")) == s("|10>")
|
||
# Tests SWAP on far-away gates
|
||
assert SWAP.on(s("|001>"), which_qbit=[0, 2]) == s("|100>")
|
||
|
||
# test Toffoli gate
|
||
assert TOFF.on(s("|000>")) == s("|000>")
|
||
assert TOFF.on(s("|100>")) == s("|100>")
|
||
assert TOFF.on(s("|110>")) == s("|111>")
|
||
assert TOFF.on(s("|111>")) == s("|110>")
|
||
|
||
# controls on 0, 1 target on 3 for 4
|
||
assert TOFF.on(s("|1100>"), which_qbit=[0, 1, 3]) == s("|1101>")
|
||
# controls on 0, 2 target on 3 for 4
|
||
assert TOFF.on(s("|1010>"), which_qbit=[0, 2, 3]) == s("|1011>")
|
||
assert TOFF.on(s("|1011>"), which_qbit=[0, 2, 3]) == s("|1010>")
|
||
assert TOFF.on(s("|1111>"), which_qbit=[0, 2, 3]) == s("|1110>")
|
||
# controls on 0, 2 target on 4 for 5
|
||
assert TOFF.on(s("|10101>"), which_qbit=[0, 2, 4]) == s("|10100>")
|
||
assert TOFF.on(s("|10111>"), which_qbit=[0, 2, 4]) == s("|10110>")
|
||
assert TOFF.on(s("|11111>"), which_qbit=[0, 2, 4]) == s("|11110>")
|
||
assert TOFF.on(s("|10100>"), which_qbit=[0, 2, 4]) == s("|10101>")
|
||
|
||
|
||
def test_sequential_measurements():
|
||
st = H.on(s("|0>"))
|
||
meas1 = st.measure()
|
||
assert st.measurement_result is not None
|
||
assert st.last_basis == computational_basis
|
||
|
||
# Measuring more times should return the same result
|
||
for i in range(10):
|
||
meas2 = st.measure()
|
||
assert meas2 == meas1
|
||
|
||
# Measuring in different basis might not yield the same result
|
||
meas3 = st.measure(basis=plus_minus_basis)
|
||
assert st.measurement_result is not None
|
||
assert st.last_basis == plus_minus_basis
|
||
|
||
# But measuring more times should still return the same result
|
||
for i in range(10):
|
||
meas4 = st.measure(basis=plus_minus_basis)
|
||
assert meas4 == meas3
|
||
|
||
|
||
def test():
|
||
# Test properties of Hilbert vector space
|
||
# The four postulates of Quantum Mechanics
|
||
# I: States | Associated to any physical system is a complex vector space
|
||
# known as the state space of the system. If the system is
|
||
# closed
|
||
# then the system is described completely by its state vector
|
||
# which is a unit vector in the space.
|
||
# Mathematically, this vector space is also a function space
|
||
assert _0 + _1 == _1 + _0 # commutativity of vector addition
|
||
assert _0 + (_1 + _p) == (_0 + _1) + _p # associativity of vector addition
|
||
assert 8 * (
|
||
_0 + _1) == 8 * _0 + 8 * _1 # Linear when multiplying by
|
||
# constants
|
||
assert _0.inner(_0) == 1 # parallel have 1 product
|
||
assert _0.inner(_1) == 0 # orthogonal have 0 product
|
||
assert _0.is_orthogonal(_1)
|
||
assert _1.inner(8 * _0) == 8 * _1.inner(
|
||
_0) # Inner product is linear multiplied by constants
|
||
assert _p.inner(_1 + _0) == _p.inner(_1) + _p.inner(
|
||
_0) # Inner product is linear in superpos of vectors
|
||
|
||
assert np.isclose(_1.length(),
|
||
1.0) # all of the vector lengths are normalized
|
||
assert np.isclose(_0.length(), 1.0)
|
||
assert np.isclose(_p.length(), 1.0)
|
||
|
||
assert _0.inner(_1) == _1.inner(
|
||
_0).complex_conjugate() # non-commutative inner product
|
||
|
||
test_to_from_angles()
|
||
|
||
# II: Dynamics | The evolution of a closed system is described by a
|
||
# unitary transformation
|
||
#
|
||
test_linear_operator()
|
||
|
||
test_eigenstuff()
|
||
|
||
# Understanding the difference between unitary and hermitians
|
||
test_unitary_hermitian()
|
||
|
||
# Pauli X gate flips the |0> to |1> and the |1> to |0>
|
||
assert X.on(_1) == _0
|
||
assert X.on(_0) == _1
|
||
|
||
# Test the Y Pauli operator with complex number literal notation
|
||
assert Y.on(_0) == State([[0],
|
||
[1j]])
|
||
assert Y.on(_1) == State([[-1j],
|
||
[0]])
|
||
|
||
# Test Pauli rotation gates
|
||
testRotPauli()
|
||
|
||
# Test 2+ qbit gates and partials
|
||
test_partials()
|
||
|
||
# III: Measurement | A quantum measurement is described by an orthonormal
|
||
# basis |e_j>
|
||
# for state space. If the initial state of the system
|
||
# is |ψ>
|
||
# then we get outcome j with probability pr(j) =
|
||
# |<e_j|ψ>|^2
|
||
# Note: The postulates are applicable on closed, isolated systems.
|
||
# Systems that are closed and are described by unitary time
|
||
# evolution by a Hamiltonian
|
||
# can be measured by projective measurements. Systems are not
|
||
# closed in reality and hence
|
||
# are immeasurable using projective measurements.
|
||
# POVM (Positive Operator-Valued Measure) is a restriction on the
|
||
# projective measurements,
|
||
# such that it encompasses everything except the environment.
|
||
|
||
assert _0.get_prob(0) == 1 # Probability for |0> in 0 is 1
|
||
assert _0.get_prob(1) == 0 # Probability for |0> in 1 is 0
|
||
|
||
assert _1.get_prob(0) == 0 # Probability for |1> in 0 is 0
|
||
assert _1.get_prob(1) == 1 # Probability for |1> in 1 is 1
|
||
|
||
assert np.isclose(_p.get_prob(0), 0.5) # Probability for |+> in 0 is 0.5
|
||
assert np.isclose(_p.get_prob(1), 0.5) # Probability for |+> in 1 is 0.5
|
||
|
||
assert _0.measure() == _0
|
||
assert _1.measure() == _1
|
||
|
||
assert s("|10>").measure() == s("|10>")
|
||
assert s("|10>").measure_partial(1) == (_1, _0)
|
||
assert s("|10>").measure_partial(2) == (_0, _1)
|
||
|
||
# measure in arbitrary basis
|
||
_0.measure(basis=[_p, _m])
|
||
s("|00>").measure(basis=[b_phi_p, b_phi_m, b_psi_m, b_psi_p])
|
||
|
||
# Maximally entangled
|
||
result, pms = b_phi_p.measure_partial(1)
|
||
if result == _0:
|
||
assert pms == _0
|
||
elif result == _1:
|
||
assert pms == _1
|
||
|
||
# Test measurement operators
|
||
test_measurement_ops()
|
||
|
||
test_sequential_measurements()
|
||
|
||
# IV: Compositing | The state space of a composite physical system
|
||
# is the tensor product of the state spaces
|
||
# of the component physical systems.
|
||
# i.e. if we have systems numbered 1 through n, ψ_i,
|
||
# then the joint state is |ψ_1> ⊗ |ψ_2> ⊗ ... ⊗ |ψ_n>
|
||
assert _0 * _0 == State([[1], [0], [0], [0]])
|
||
assert _0 * _1 == State([[0], [1], [0], [0]])
|
||
assert _1 * _0 == State([[0], [0], [1], [0]])
|
||
assert _1 * _1 == State([[0], [0], [0], [1]])
|
||
|
||
# CNOT applies a control qubit on a target.
|
||
# If the control is a |0>, target remains unchanged.
|
||
# If the control is a |1>, target is flipped.
|
||
assert CNOT.on(_0 * _0) == _0 * _0
|
||
assert CNOT.on(_0 * _1) == _0 * _1
|
||
assert CNOT.on(_1 * _0) == _1 * _1
|
||
assert CNOT.on(_1 * _1) == _1 * _0
|
||
|
||
# ALL FOUR NOW - Create a Bell state (I) that has H|0> (II), measures (
|
||
# III) and composition in CNOT (IV)
|
||
# Bell state
|
||
# First - create a superposition
|
||
superpos = H.on(_0)
|
||
assert superpos == _p
|
||
|
||
# Then CNOT the superposition with a |0> qubit
|
||
bell = CNOT.on(superpos * _0)
|
||
assert bell == State([[1 / np.sqrt(2)],
|
||
[0.],
|
||
[0.],
|
||
[1 / np.sqrt(2)], ])
|
||
|
||
assert np.isclose(bell.get_prob(0b00),
|
||
0.5) # Probability for bell in 00 is 0.5
|
||
assert np.isclose(bell.get_prob(0b01),
|
||
0) # Probability for bell in 01 is 0.
|
||
assert np.isclose(bell.get_prob(0b10),
|
||
0) # Probability for bell in 10 is 0.
|
||
assert np.isclose(bell.get_prob(0b11),
|
||
0.5) # Probability for bell in 11 is 0.5
|
||
|
||
################################
|
||
assert _0.x(_0) == Matrix([[1, 0],
|
||
[0, 0]])
|
||
assert _0.x(_1) == Matrix([[0, 1],
|
||
[0, 0]])
|
||
assert _1.x(_0) == Matrix([[0, 0],
|
||
[1, 0]])
|
||
assert _1.x(_1) == Matrix([[0, 0],
|
||
[0, 1]])
|
||
|
||
|
||
def test_to_from_angles():
|
||
for q in [_0, _1, _p, _m]:
|
||
angles = q.to_bloch_angles()
|
||
s = State.from_bloch_angles(*angles)
|
||
assert q == s
|
||
|
||
assert State.from_bloch_angles(0, 0) == _0
|
||
assert State.from_bloch_angles(np.pi, 0) == _1
|
||
assert State.from_bloch_angles(np.pi / 2, 0) == _p
|
||
assert State.from_bloch_angles(np.pi / 2, np.pi) == _m
|
||
for theta, phi, qbit in [
|
||
(0, 0, _0),
|
||
(np.pi, 0, _1),
|
||
(np.pi / 2, 0, _p),
|
||
(np.pi / 2, np.pi, _m),
|
||
]:
|
||
s = State.from_bloch_angles(theta, phi)
|
||
assert s == qbit
|
||
|
||
|
||
def naive_load_test(N):
|
||
import os
|
||
import psutil
|
||
import gc
|
||
from time import time
|
||
from sys import getsizeof
|
||
|
||
print(
|
||
"{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format(
|
||
"qbits",
|
||
"kron_len",
|
||
"mem_used",
|
||
"mem_per_q",
|
||
"getsizeof",
|
||
"getsiz/len",
|
||
"nbytes",
|
||
"nbytes/len",
|
||
"time"))
|
||
|
||
_0 = State([[1], [0]], name='0')
|
||
|
||
process = psutil.Process(os.getpid())
|
||
mem_init = process.memory_info().rss
|
||
for i in range(2, N + 1):
|
||
start = time()
|
||
m = _0
|
||
for _ in range(i):
|
||
m = m * _0
|
||
len_m = len(m)
|
||
elapsed = time() - start
|
||
mem_b = process.memory_info().rss - mem_init
|
||
print(
|
||
"{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {"
|
||
":>10}".format(
|
||
i,
|
||
len(m),
|
||
sizeof_fmt(mem_b),
|
||
sizeof_fmt(mem_b / len_m),
|
||
sizeof_fmt(getsizeof(m)),
|
||
sizeof_fmt(getsizeof(m) / len_m),
|
||
sizeof_fmt(m.m.nbytes),
|
||
sizeof_fmt(m.m.nbytes / len_m),
|
||
np.round(elapsed, 2)))
|
||
|
||
gc.collect()
|
||
|
||
|
||
class QuantumCircuit(object):
|
||
def __init__(self, n_qubits: int, initial_steps=None):
|
||
self.n_qubits = n_qubits
|
||
if not initial_steps:
|
||
self.steps = [[_0 for _ in range(n_qubits)], ]
|
||
else:
|
||
self.steps = initial_steps
|
||
self._called_add_row = 0
|
||
|
||
@property
|
||
def rows(self):
|
||
return self._called_add_row
|
||
|
||
def add_row_step(self, row: int, step: int, qbit_state):
|
||
if len(self.steps) <= step:
|
||
self.steps += [[I for _ in range(self.n_qubits)] for _ in
|
||
range(len(self.steps) - step + 1)]
|
||
self.steps[step][row] = qbit_state
|
||
|
||
def add_step(self, step_data: list):
|
||
if len(step_data) != self.n_qubits:
|
||
raise RuntimeError(
|
||
"Length of step is: {}, should be: {}".format(len(step_data),
|
||
self.n_qubits))
|
||
step_i = len(step_data)
|
||
for row, qubit_state in enumerate(step_data):
|
||
self.add_row_step(row, step_i, qubit_state)
|
||
|
||
def add_steps(self, steps_data: list):
|
||
for step_data in steps_data:
|
||
self.add_step(step_data)
|
||
|
||
def add_row(self, row_data: list):
|
||
if self._called_add_row >= self.n_qubits:
|
||
raise RuntimeError("Adding more rows than qubits")
|
||
for step_i, qubit_state in enumerate(row_data):
|
||
self.add_row_step(self._called_add_row, step_i + 1, qubit_state)
|
||
self._called_add_row += 1
|
||
|
||
def add_rows(self, rows_data: list):
|
||
for row_data in rows_data:
|
||
self.add_row(row_data)
|
||
|
||
def compose_quantum_state(self, step):
|
||
return reduce((lambda x, y: x * y), step)
|
||
|
||
def step(self):
|
||
if self.current_step == 0 and self.current_state == self.HALT_STATE:
|
||
self.current_state = self.RUNNING_STATE
|
||
if self.current_step >= len(self.steps):
|
||
self.current_state = self.HALT_STATE
|
||
raise RuntimeWarning("Halted")
|
||
running_step = self.steps[self.current_step]
|
||
step_quantum_state = self.compose_quantum_state(running_step)
|
||
if not self.current_quantum_state:
|
||
self.current_quantum_state = step_quantum_state
|
||
else:
|
||
self.current_quantum_state = step_quantum_state.inner(
|
||
self.current_quantum_state)
|
||
self.current_step += 1
|
||
|
||
def run(self):
|
||
for _ in self.steps:
|
||
self.step()
|
||
self.current_state = self.HALT_STATE
|
||
|
||
def reset(self):
|
||
self.current_step = 0
|
||
self.current_quantum_state = None
|
||
self.current_state = self.HALT_STATE
|
||
|
||
def print(self):
|
||
print("=" * 3 * len(self.steps))
|
||
for line_no in range(self.n_qubits):
|
||
line = ''
|
||
for step in self.steps:
|
||
state_repr = repr(step[line_no])
|
||
line += state_repr
|
||
print(line)
|
||
print("=" * 3 * len(self.steps))
|
||
|
||
|
||
class QuantumProcessor(object):
|
||
HALT_STATE = "HALT"
|
||
RUNNING_STATE = "RUNNING"
|
||
|
||
def __init__(self, circuit: QuantumCircuit):
|
||
if circuit.rows != circuit.n_qubits:
|
||
raise Exception(
|
||
"Declared circuit with n_qubits: {} but called add_row: {"
|
||
"}".format(
|
||
circuit.n_qubits, circuit.rows
|
||
))
|
||
self.circuit = circuit
|
||
self.c_step = 0
|
||
self.c_q_state = None
|
||
self.c_state = self.HALT_STATE
|
||
self.reset()
|
||
|
||
def compose_quantum_state(self, step):
|
||
return reduce((lambda x, y: x * y), step)
|
||
|
||
def step(self):
|
||
if self.c_step == 0 and self.c_state == self.HALT_STATE:
|
||
self.c_state = self.RUNNING_STATE
|
||
if self.c_step >= len(self.circuit.steps):
|
||
self.c_state = self.HALT_STATE
|
||
raise RuntimeWarning("Halted")
|
||
step_quantum_state = self.get_next_step()
|
||
if not self.c_q_state:
|
||
self.c_q_state = step_quantum_state
|
||
else:
|
||
self.c_q_state = State((step_quantum_state.inner(self.c_q_state)).m)
|
||
self.c_step += 1
|
||
|
||
def get_next_step(self):
|
||
running_step = self.circuit.steps[self.c_step]
|
||
step_quantum_state = self.compose_quantum_state(running_step)
|
||
return step_quantum_state
|
||
|
||
def run(self):
|
||
for _ in self.circuit.steps:
|
||
self.step()
|
||
self.c_state = self.HALT_STATE
|
||
|
||
def reset(self):
|
||
self.c_step = 0
|
||
self.c_q_state = None
|
||
self.c_state = self.HALT_STATE
|
||
|
||
def measure(self):
|
||
if self.c_state != self.HALT_STATE:
|
||
raise RuntimeError("Processor is still running")
|
||
return self.c_q_state.measure()
|
||
|
||
def run_n(self, n: int):
|
||
for i in range(n):
|
||
self.run()
|
||
result = self.measure()
|
||
print("Run {}: {}".format(i, result))
|
||
self.reset()
|
||
|
||
def get_sample(self, n: int):
|
||
rv = defaultdict(int)
|
||
for i in range(n):
|
||
self.run()
|
||
result = self.measure()
|
||
rv[result] += 1
|
||
self.reset()
|
||
return rv
|
||
|
||
@staticmethod
|
||
def print_sample(rv):
|
||
for k, v in rv.items():
|
||
print("{}: {}".format(k, v))
|
||
|
||
|
||
def test_quantum_processor():
|
||
# Produce Bell state between 0 and 2 qubit
|
||
qc = QuantumCircuit(3)
|
||
qc.add_row([H, C_partial])
|
||
qc.add_row([_, _])
|
||
qc.add_row([_, x_partial])
|
||
qc.print()
|
||
qp = QuantumProcessor(qc)
|
||
qp.print_sample(qp.get_sample(100))
|
||
|
||
|
||
def test_light():
|
||
# http://alienryderflex.com/polarizer/
|
||
hor_filter = MeasurementOperator.create_from_basis(Matrix(_0.m), name='h')
|
||
diag_filter = MeasurementOperator.create_from_basis(Matrix(_p.m), name='d')
|
||
ver_filter = MeasurementOperator.create_from_basis(Matrix(_1.m), name='v')
|
||
|
||
def random_light():
|
||
random_pol = Vector(
|
||
[[np.random.uniform(0, 1)], [np.random.uniform(0, 1)]])
|
||
return State.normalize(random_pol)
|
||
|
||
def experiment(id, random_ls, filters):
|
||
total_light = defaultdict(int)
|
||
human_state = "Light"
|
||
for filter in filters:
|
||
human_state += "->{}".format(filter.name)
|
||
|
||
for r in random_ls:
|
||
st = filters[0].on(r)
|
||
for filter in filters:
|
||
st = filter.on(st)
|
||
st = State(st, allow_unnormalized=True)
|
||
total_light[st.measure(allow_empty=True)] += 1
|
||
print("{}. {}:\n {}".format(id, human_state, dict(total_light)))
|
||
|
||
its = 100
|
||
random_lights = [random_light() for _ in range(its)]
|
||
|
||
# Just horizonal - should result in some light
|
||
experiment(1, random_lights, [hor_filter, ])
|
||
|
||
# Vertical after horizonal - should result in 0
|
||
experiment(2, random_lights, [ver_filter, hor_filter])
|
||
|
||
# TODO: Something is wrong here...
|
||
# Vertical after diagonal after horizontal - should result in ~ 50%
|
||
# compared to only Horizontal
|
||
experiment(3, random_lights, [ver_filter, diag_filter, hor_filter])
|
||
|
||
|
||
def get_bin_fmt(count):
|
||
return "{:0" + str(int(np.ceil(np.log2(count)))) + "b}"
|
||
|
||
|
||
def generate_bin(i, count):
|
||
format_str = get_bin_fmt(count)
|
||
return format_str.format(i)
|
||
|
||
|
||
def generate_bins(count):
|
||
format_str = get_bin_fmt(count)
|
||
return [format_str.format(i) for i in range(count)]
|
||
|
||
|
||
if __name__ == "__main__":
|
||
test()
|
||
test_quantum_processor()
|
||
test_light()
|