diff --git a/lib_q_computer_math.py b/lib_q_computer_math.py index 5afee9c..474704c 100644 --- a/lib_q_computer_math.py +++ b/lib_q_computer_math.py @@ -8,13 +8,26 @@ 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] +REPR_TENSOR_OP = "⊗" +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): - """Represents a quantum state as a vector in Hilbert space""" + """Wraps a Matrix... it's for my understanding, this could easily probably be np.array""" def __init__(self, m: ListOrNdarray = None, *args, **kwargs): """ @@ -32,12 +45,12 @@ 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): return bool(np.allclose(self.m, other)) - elif isinstance(other, TwoQubitGate): + elif isinstance(other, TwoQubitPartial): return False return np.allclose(self.m, other.m) @@ -47,9 +60,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 +74,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 +102,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, *args, **kwargs) + 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""" - super().__init__(m) + """State vector representing quantum state""" + super().__init__(m, *args, **kwargs) 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) @classmethod def from_angles(cls, theta, phi): @@ -125,7 +161,16 @@ class State(Matrix): def __repr__(self): if self.name: return '|{}>'.format(self.name) - return str(self.m) + for well_known_state in well_known_states: + if self == well_known_state: + return repr(well_known_state) + for used_state in UNIVERSE_STATES: + if self.m == used_state: + return repr(used_state) + next_state_sub = ''.join([REPR_MATH_SUBSCRIPT_NUMBERS[int(d)] for d in str(len(UNIVERSE_STATES))]) + state_name = '|{}{}>'.format(REPR_GREEK_PSI, next_state_sub) + UNIVERSE_STATES.append(self.m) + return state_name def norm(self): """Norm/Length of the vector = sqrt()""" @@ -164,7 +209,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) @@ -172,22 +250,22 @@ 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): + return State(np.dot(self.m, other.m)) def __repr__(self): if self.name: @@ -195,6 +273,68 @@ class UnitaryOperator(Operator, UnitaryMatrix): return str(self.m) +# 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) +# +# 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 TwoQubitPartial(object): + def __init__(self, rpr): + self.rpr = rpr + self.operator = None + + def __repr__(self): + return str("-{}-".format(self.rpr)) + + +C_ = TwoQubitPartial("C") +x_ = TwoQubitPartial("x") + + +class TwoQubitOperator(UnitaryOperator): + def __init__(self, m: ListOrNdarray, A: TwoQubitPartial, B: TwoQubitPartial, + 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) + + """ Define States and Operators """ @@ -207,9 +347,13 @@ _1 = State([[0], [1]], name='1') -_p = State([[1 / np.sqrt(2)], [1 / np.sqrt(2)]], name='+') +_p = State([[1 / np.sqrt(2)], + [1 / np.sqrt(2)]], + name='+') -_m = State([[1 / np.sqrt(2)], [-1 / np.sqrt(2)]], name='-') +_m = State([[1 / np.sqrt(2)], + [- 1 / np.sqrt(2)]], + name='-') _00 = State([[1], [0], @@ -223,6 +367,8 @@ _11 = State([[0], [1]], name='11') +well_known_states = [_0, _1, _p, _m] + _ = I = UnitaryOperator([[1, 0], [0, 1]], name="-") @@ -243,39 +389,16 @@ 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], ]) +CNOT = TwoQubitOperator([[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + [0, 0, 1, 0], ], + TwoQubitPartial("C"), + TwoQubitPartial("x"), + I, + X) - -# 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) -# -# 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 TwoQubitGate(object): - def __init__(self, rpr): - self.rpr = rpr - - def __repr__(self): - return str("-{}-".format(self.rpr)) - - -C = TwoQubitGate("C") -x = TwoQubitGate("x") +C, x = CNOT.A, CNOT.B # TODO: End Hacky way to define 2-qbit gate @@ -283,8 +406,6 @@ x = TwoQubitGate("x") def test(): - _p = State([[1 / np.sqrt(2)], [1 / np.sqrt(2)]], name='+') - # Test properties of Hilbert vector space # The four postulates of Quantum Mechanics # I: States | Associated to any physical system is a complex vector space @@ -448,18 +569,11 @@ def naive_load_test(N): gc.collect() -class QuantumProcessor(object): - HALT_STATE = "HALT" - RUNNING_STATE = "RUNNING" - +class QuantumCircuit(object): def __init__(self, n_qubits: int): self.n_qubits = n_qubits self.steps = [[_0 for _ in range(n_qubits)], ] self._called_add_row = 0 - self.current_step = 0 - self.current_quantum_state = None - self.current_state = self.HALT_STATE - self.reset() def add_row_step(self, row: int, step: int, qbit_state): if len(self.steps) <= step: @@ -543,10 +657,60 @@ class QuantumProcessor(object): print(line) print("=" * 3 * len(self.steps)) + +class QuantumProcessor(object): + HALT_STATE = "HALT" + RUNNING_STATE = "RUNNING" + + def __init__(self, circuit: QuantumCircuit): + 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): + if any([type(s) is TwoQubitPartial for s in step]): + two_qubit_gates = filter(lambda s: type(s) is TwoQubitPartial, 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): + 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 | 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.current_state != self.HALT_STATE: + if self.c_state != self.HALT_STATE: raise RuntimeError("Processor is still running") - return self.current_quantum_state.measure() + return self.c_q_state.measure() def run_n(self, n: int): for i in range(n): @@ -569,19 +733,19 @@ class QuantumProcessor(object): def test_quantum_processor(): # Produce Bell state between 0 and 2 qubit - qp = QuantumProcessor(3) - qp.add_row([H, C]) - qp.add_row([_, _]) - qp.add_row([_, x]) - qp.print() + qc = QuantumCircuit(3) + qc.add_row([H, C]) + qc.add_row([_, _]) + qc.add_row([_, x]) + qc.print() + qp = QuantumProcessor(qc) qp.get_sample(100) if __name__ == "__main__": test() test_to_from_angles() - - # test_quantum_processor() + test_quantum_processor() # print(_1 | _0) # qubit = _00 + _11 # print(qubit)