diff --git a/lib_q_computer_math.py b/lib_q_computer_math.py index 2cc657b..adb546f 100644 --- a/lib_q_computer_math.py +++ b/lib_q_computer_math.py @@ -8,6 +8,8 @@ from typing import Union import numpy as np from numbers import Complex +import sympy + from load_test import sizeof_fmt ListOrNdarray = Union[list, np.ndarray] @@ -32,7 +34,7 @@ class Matrix(object): raise TypeError("m needs to be a list or ndarray type") def __add__(self, other): - return other.__class__(self.m + other.m) + return Matrix(self.m + other.m) def __eq__(self, other): if isinstance(other, Complex): @@ -47,9 +49,9 @@ class Matrix(object): with another state is a composition via the Kronecker/Tensor product """ if isinstance(other, Complex): - return self.__class__(self.m * other) + return Matrix(self.m * other) elif isinstance(other, Matrix): - return other.__class__(np.kron(self.m, other.m)) + return self.__class__(np.kron(self.m, other.m)) raise NotImplementedError def __rmul__(self, other): @@ -61,14 +63,22 @@ class Matrix(object): def __or__(self, other): """Define inner product: """ - return other.__class__(np.dot(self._conjugate_transpose(), other.m)) + m = np.dot(self._conjugate_transpose(), other.m) + try: + return self.__class__(m) + except: + try: + return other.__class__(m) + except: + ... + return Matrix(m) def __len__(self): return len(self.m) def outer(self, other): """Define outer product |0><0|""" - return self.__class__(np.outer(self.m, other.m)) + 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)""" @@ -81,17 +91,32 @@ class Matrix(object): return self.m.conjugate() def conjugate_transpose(self): - return self.__class__(self._conjugate_transpose()) + return Matrix(self._conjugate_transpose()) def complex_conjugate(self): - return self.__class__(self._complex_conjugate()) + return Matrix(self._complex_conjugate()) -class State(Matrix): +class Vector(Matrix): + def __init__(self, m: ListOrNdarray = None, *args, **kwargs): + super().__init__(m) + if not self._is_vector(): + raise TypeError("Not a vector") + + def _is_vector(self): + return self.m.shape[1] == 1 + + +class State(Vector): def __init__(self, m: ListOrNdarray = None, name: str = '', *args, **kwargs): - """An Operator turns one function into another""" + """State vector representing quantum state""" super().__init__(m) self.name = name + if not self._is_normalized(): + raise TypeError("Not a normalized state vector") + + def _is_normalized(self): + return np.isclose(np.sum(np.abs(self.m ** 2)), 1.0) def __repr__(self): if self.name: @@ -135,7 +160,40 @@ class Operator(object): return self.func(*args, **kwargs) -class UnitaryMatrix(Matrix): +class LinearOperator(Operator): + def __init__(self, func=None, *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): + # TODO: How to verify if the func is linear? + # in case of Unitary Operator, self.func is a lambda that takes a Matrix (assumes has .m component) + return True + # a, b = sympy.symbols('a, b') + # expr, vars_ = a+b, [a, b] + # for x in vars_: + # for y in vars_: + # try: + # if not sympy.Eq(sympy.diff(expr, x, y), 0): + # return False + # except TypeError: + # return False + # return True + + +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 UnitaryMatrix(SquareMatrix): def __init__(self, m: ListOrNdarray, *args, **kwargs): """Represents a Unitary matrix that satisfies UU+ = I""" super().__init__(m, *args, **kwargs) @@ -143,22 +201,26 @@ class UnitaryMatrix(Matrix): raise TypeError("Not a Unitary matrix") def _is_unitary(self): - """Checks if the matrix is - 1. square and - 2. the product of itself with conjugate transpose is Identity UU+ = I""" - is_square = self.m.shape[0] == self.m.shape[1] + """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 is_square and np.isclose(UU_, I).all() + return np.isclose(UU_, I).all() -class UnitaryOperator(Operator, UnitaryMatrix): +class UnitaryOperator(LinearOperator, UnitaryMatrix): def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs): - """UnitaryOperator inherits from both Operator and a Unitary matrix + """UnitaryOperator inherits from both LinearOperator and a Unitary matrix It is used to act on a State vector by defining the operator to be the dot product""" self.name = name UnitaryMatrix.__init__(self, m=m, *args, **kwargs) - Operator.__init__(self, func=lambda other: State(np.dot(self.m, other.m)), *args, **kwargs) + LinearOperator.__init__(self, func=self.operator_func, *args, **kwargs) + + def operator_func(self, other): + if not hasattr(other, "m"): + other_m = other + else: + other_m = other.m + return State(np.dot(self.m, other_m)) def __repr__(self): if self.name: @@ -198,11 +260,6 @@ H = UnitaryOperator([[1 / np.sqrt(2), 1 / np.sqrt(2)], [1 / np.sqrt(2), -1 / np.sqrt(2)], ], name="H") -CNOT = UnitaryOperator([[1, 0, 0, 0], - [0, 1, 0, 0], - [0, 0, 0, 1], - [0, 0, 1, 0], ]) - # TODO - 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) @@ -224,13 +281,59 @@ CNOT = UnitaryOperator([[1, 0, 0, 0], class TwoQubitGate(object): def __init__(self, rpr): self.rpr = rpr + self.operator = None def __repr__(self): return str("-{}-".format(self.rpr)) -C = TwoQubitGate("C") -x = TwoQubitGate("x") +C_ = TwoQubitGate("C") +x_ = TwoQubitGate("x") + + +class TwoQubitOperator(UnitaryOperator): + def __init__(self, m: ListOrNdarray, A: TwoQubitGate, B: TwoQubitGate, + 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 + + def verify_step(self, step): + if not (step.count(self.A) == 1 and step.count(self.B) == 1): + raise RuntimeError("Both CONTROL and TARGET need to be defined in the same step exactly once") + + def compose(self, step, state): + # TODO: Hacky way to do CNOT + # Should generalize for a 2-Qubit gate + # _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) + + +CNOT = TwoQubitOperator([[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + [0, 0, 1, 0], ], + TwoQubitGate("C"), + TwoQubitGate("x"), + I, + X) + +C, x = CNOT.A, CNOT.B + # TODO: End Hacky way to define 2-qbit gate ########################################################### @@ -423,24 +526,13 @@ class QuantumProcessor(object): self.add_row(row_data) def compose_quantum_state(self, step): - if C in step or x in step: - if not (step.count(C) == 1 and step.count(x) == 1): - raise RuntimeError("Both CONTROL and TARGET need to be defined in the same step exactly once") - # TODO: Hacky way to do CNOT - # Should generalize for a 2-Qubit gate - # _0.x(_0) * Matrix(I.m) + _1.x(_1) * Matrix(X.m) - outer_0, outer_1 = [], [] - for s in step: - if s == C: - outer_0.append(_0.x(_0)) - outer_1.append(_1.x(_1)) - elif s == x: - outer_0.append(Matrix(I.m)) - outer_1.append(Matrix(X.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) + if any([type(s) is TwoQubitGate for s in step]): + two_qubit_gates = filter(lambda s: type(s) is TwoQubitGate, 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): @@ -454,7 +546,7 @@ class QuantumProcessor(object): if not self.current_quantum_state: self.current_quantum_state = step_quantum_state else: - self.current_quantum_state = step_quantum_state | self.current_quantum_state + self.current_quantum_state = State((step_quantum_state | self.current_quantum_state).m) self.current_step += 1 def run(self):