From 60a11f143ef6580b144a31c13f316ab7bcea0f02 Mon Sep 17 00:00:00 2001 From: Daniel Tsvetkov Date: Sat, 28 Mar 2020 22:17:40 +0100 Subject: [PATCH] partial application of multiqubit gates --- lib.py | 505 ++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 301 insertions(+), 204 deletions(-) diff --git a/lib.py b/lib.py index 875d333..3af7e11 100644 --- a/lib.py +++ b/lib.py @@ -13,6 +13,7 @@ ListOrNdarray = Union[list, np.ndarray] REPR_EMPTY_SET = "Ø" REPR_TENSOR_OP = "⊗" +REPR_TARGET = "⊕" REPR_GREEK_BETA = "β" REPR_GREEK_PSI = "ψ" REPR_GREEK_PHI = "φ" @@ -26,7 +27,8 @@ UNIVERSE_STATES = [] class Matrix(object): - """Wraps a Matrix... it's for my understanding, this could easily probably be np.array""" + """Wraps a Matrix... it's for my understanding, this could easily + probably be np.array""" def __init__(self, m: ListOrNdarray = None, *args, **kwargs): """ @@ -52,8 +54,6 @@ class Matrix(object): def __eq__(self, other): if isinstance(other, Complex): return bool(np.allclose(self.m, other)) - elif isinstance(other, PartialQubit): - return False return np.allclose(self.m, other.m) def __mul_or_kron__(self, other): @@ -182,7 +182,9 @@ class State(Vector): 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)") + 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 @@ -193,7 +195,8 @@ class State(Vector): div = np.sin(theta / 2) if div == 0: - # here is doesn't matter what phi is as phase at the poles is arbitrary + # here is doesn't matter what phi is as phase at the poles is + # arbitrary phi = 0 else: exp = m1 / div @@ -234,10 +237,13 @@ class State(Vector): return repr(used_state) except: ... - next_state_sub = ''.join([REPR_MATH_SUBSCRIPT_NUMBERS[int(d)] for d in str(len(UNIVERSE_STATES))]) + 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) - # matrix_rep = "{}".format(self.m).replace('[', '').replace(']', '').replace('\n', '|').strip() + # matrix_rep = "{}".format(self.m).replace('[', '').replace(']', + # '').replace('\n', '|').strip() # state_name = '|{}> = {}'.format(self.name, matrix_rep) state_name = '|{}>'.format(self.name) return state_name @@ -251,7 +257,8 @@ class State(Vector): 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""" + """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): @@ -284,8 +291,10 @@ class State(Vector): 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? + 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") """ # if self.measurement_result: @@ -320,7 +329,8 @@ class State(Vector): weights = list(np.array(weights) / sum(weights)) format_str = self.get_fmt_of_element() - choices = empty_choices + [format_str.format(i) for i in range(len(weights))] + choices = empty_choices + [format_str.format(i) for i in + range(len(weights))] weights = empty_weights + weights self.measurement_result = random.choices(choices, weights)[0] @@ -329,7 +339,8 @@ class State(Vector): 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 + 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) @@ -341,7 +352,10 @@ class State(Vector): """ 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)) + 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'] @@ -352,12 +366,15 @@ class State(Vector): # [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] + 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]])) + normalization_factor = np.sqrt( + sum([np.abs(i) ** 2 for i in weights[measurement_result]])) new_m = weights[measurement_result] / normalization_factor return str(measurement_result), State(new_m.reshape((len(new_m), 1))) @@ -375,7 +392,9 @@ class State(Vector): 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") + raise Exception( + "State from string should be of the form |00..01> with all " + "numbers either 0 or 1") return cls(arr) @staticmethod @@ -405,7 +424,8 @@ def test_measure_partial(): def normalize_state(vector: Vector): - """Normalize a state by dividing by the square root of sum of the squares of states""" + """Normalize a state by dividing by the square root of sum of the squares + of states""" norm_coef = np.sqrt(np.sum(np.array(vector.m) ** 2)) if norm_coef == 0: raise TypeError("zero-sum vector") @@ -464,11 +484,11 @@ class Operator(object): """An Operator turns one function into another""" self.func = func - def on(self, *args): - return self(*args) + def on(self, *args, **kwargs): + return self(*args, **kwargs) - def __call__(self, *args): - return self.func(*args) + def __call__(self, *args, **kwargs): + return self.func(*args, **kwargs) class LinearOperator(Operator): @@ -480,8 +500,10 @@ class LinearOperator(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 + # 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_: @@ -501,7 +523,8 @@ def test_linear_operator(): assert _times_2.on(5) == 10 assert _times_2(5) == 10 - assert_raises(TypeError, "Not a linear operator", LinearOperator, sp.Lambda(_x, _x ** 2)) + assert_raises(TypeError, "Not a linear operator", LinearOperator, + sp.Lambda(_x, _x ** 2)) class SquareMatrix(Matrix): @@ -515,9 +538,12 @@ class SquareMatrix(Matrix): 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""" + 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 @@ -525,7 +551,8 @@ class LinearTransformation(LinearOperator, Matrix): def _is_linear(self): # Every matrix transformation is a linear transformation - # https://www.mathbootcamps.com/proof-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): @@ -555,7 +582,8 @@ class UnitaryMatrix(SquareMatrix): raise TypeError("Not a Unitary matrix") def _is_unitary(self): - """Checks if the matrix product of itself with conjugate transpose is Identity UU+ = I""" + """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() @@ -574,25 +602,67 @@ class HermitianMatrix(SquareMatrix): class UnitaryOperator(LinearTransformation, UnitaryMatrix): - def __init__(self, m: ListOrNdarray, name: str = '', *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""" + def __init__(self, m: ListOrNdarray, name: str = '', partials=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: + raise Exception("Please define partials in a non-single operator") UnitaryMatrix.__init__(self, m=m, *args, **kwargs) - LinearTransformation.__init__(self, m=m, func=self.operator_func, *args, **kwargs) + LinearTransformation.__init__(self, m=m, func=self.operator_func, *args, + **kwargs) self.name = name + self.partials = partials 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: - 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)) + if this_cols == other_rows: + 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 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: + # 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]) + extended_m, next_partial = [[], []], 0 + for qbit in range(total_qbits): + if qbit not in which_qbit: + extended_m[0].append(I) + extended_m[1].append(I) + else: + this_partial = self.partials[next_partial] + if this_partial == C_partial: + extended_m[0].append(s("|0><0|")) + extended_m[1].append(s("|1><1|")) + else: + extended_m[0].append(I) + extended_m[1].append(this_partial) + next_partial += 1 + 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") - extended_m = [self if i == which_qbit else I for i in range(int(np.log2(other_rows)))] - extended_op = UnitaryOperator(np.prod(extended_m).m) - return State(np.dot(extended_op.m, other.m)) - return State(np.dot(self.m, other.m)) + extended_op = UnitaryOperator(new_m, name=self.name, partials=self.partials) + return State(np.dot(extended_op.m, other.m)) def __repr__(self): if self.name: @@ -608,14 +678,17 @@ class Gate(UnitaryOperator): 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""" + """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""" + """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): @@ -643,71 +716,6 @@ class DensityMatrix(HermitianOperator): return np.isclose(np.trace(self.m), 1.0) -# 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) - -class PartialQubit(object): - def __init__(self, rpr): - self.rpr = rpr - self.operator = None - - def __repr__(self): - return str("-{}-".format(self.rpr)) - - -C_partial = PartialQubit("C") -x_partial = PartialQubit("x") - - -class TwoQubitOperator(UnitaryOperator): - def __init__(self, m: ListOrNdarray, A: PartialQubit, B: PartialQubit, - A_p: UnitaryOperator, B_p: UnitaryOperator, *args, **kwargs): - super().__init__(m, *args, **kwargs) - A.operator, B.operator = self, self - self.A = A - self.B = B - self.A_p = A_p - self.B_p = B_p - - A.operator = self - B.operator = self - - def verify_step(self, step): - if not (step.count(self.A) == 1 and step.count(self.B) == 1): - raise RuntimeError("Both {} and {} need to be defined in the same step exactly once".format( - self.A, self.B - )) - - def compose(self, step, state): - # _0.x(_0) * Matrix(I.m) + _1.x(_1) * Matrix(X.m) - outer_0, outer_1 = [], [] - for s in step: - if s == self.A: - outer_0.append(_0.x(_0)) - outer_1.append(_1.x(_1)) - elif s == self.B: - outer_0.append(Matrix(self.A_p.m)) - outer_1.append(Matrix(self.B_p.m)) - else: - outer_0.append(Matrix(s.m)) - outer_1.append(Matrix(s.m)) - return reduce((lambda x, y: x * y), outer_0) + reduce((lambda x, y: x * y), outer_1) - - """ Define States and Operators """ @@ -767,73 +775,92 @@ 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="-") + [0, 1]], + name="-") X = Gate([[0, 1], - [1, 0]], - name="X") + [1, 0]], + name="X") Y = Gate([[0, -1j], - [1j, 0]], - name="Y") + [1j, 0]], + name="Y") Z = Gate([[1, 0], - [0, -1]], - name="Z") - + [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 +# 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 +# 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)], ], + [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") - + [-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") - + [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") - + [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)]], + [0, np.power(np.e, 1j * phi / 2)]], name="T") -CNOT = TwoQubitOperator([ +# 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) + +C_partial = Gate(I.m, name="C") +x_partial = Gate(X.m, name=REPR_TARGET) + +CNOT = Gate([ [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], -], C_partial, x_partial, I, X) +], name="CNOT", partials=[C_partial, x_partial]) - -# TOFFOLLI_GATE = ThreeQubitOperator([ -# [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], -# ], C_partial, C_partial, x_partial, I, I, X) +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): @@ -857,12 +884,15 @@ 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 + # 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_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 = [ @@ -876,14 +906,16 @@ def test_unitary_hermitian(): [0, 1], [1, 0], ] - assert_not_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix, u_and_h) + 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 Hermitian matrix", HermitianMatrix, + not_u_not_h) assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, not_u_not_h) @@ -947,32 +979,69 @@ def test_eigenstuff(): [(1.0, HorizontalVector([1., 0.])), (0., HorizontalVector([0., 1.]))] +def test_partials(): + # 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>")) + # 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 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("|1100>"), which_qbit=[0, 1, 3]) == s("|1101>") + + 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 + # 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 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 _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(_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 + 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 + # II: Dynamics | The evolution of a closed system is described by a + # unitary transformation # test_linear_operator() @@ -994,14 +1063,23 @@ def test(): # Test Pauli rotation gates testRotPauli() - # 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) = ||^2 + # 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) = + # ||^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 + # 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, + # 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 @@ -1052,7 +1130,8 @@ def test(): 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) + # 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) @@ -1065,10 +1144,14 @@ def test(): [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 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], @@ -1108,16 +1191,17 @@ def naive_load_test(N): 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")) + 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') @@ -1131,16 +1215,18 @@ def naive_load_test(N): 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))) + 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() @@ -1160,12 +1246,15 @@ class QuantumCircuit(object): 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 += [[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)) + 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) @@ -1187,7 +1276,8 @@ class QuantumCircuit(object): def compose_quantum_state(self, step): partials = [op for op in step if isinstance(op, PartialQubit)] - # TODO: No more than 1 TwoQubitGate **OR** UnitaryOperator can be used in a step + # TODO: No more than 1 TwoQubitGate **OR** UnitaryOperator can be + # used in a step for partial in partials: two_qubit_op = partial.operator two_qubit_op.verify_step() @@ -1205,7 +1295,8 @@ class QuantumCircuit(object): 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_quantum_state = step_quantum_state.inner( + self.current_quantum_state) self.current_step += 1 def run(self): @@ -1235,9 +1326,11 @@ class QuantumProcessor(object): 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 - )) + 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 @@ -1245,13 +1338,6 @@ class QuantumProcessor(object): self.reset() def compose_quantum_state(self, step): - if any([type(s) is PartialQubit for s in step]): - two_qubit_gates = filter(lambda s: type(s) is PartialQubit, step) - state = [] - for two_qubit_gate in two_qubit_gates: - two_qubit_gate.operator.verify_step(step) - state = two_qubit_gate.operator.compose(step, state) - return state return reduce((lambda x, y: x * y), step) def step(self): @@ -1327,7 +1413,8 @@ def test_light(): 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)]]) + random_pol = Vector( + [[np.random.uniform(0, 1)], [np.random.uniform(0, 1)]]) return normalize_state(random_pol) def experiment(id, random_ls, filters): @@ -1354,12 +1441,22 @@ def test_light(): 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 + # 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 = "{:0" + str(int(np.ceil(np.log2(count)))) + "b}" + format_str = get_bin_fmt(count) return [format_str.format(i) for i in range(count)]