Merge branch 'master' of gitlab.com:pisquared/quantum

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: <self|other>
         """
-        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"""
+        """State vector representing quantum state"""
-        super().__init__(m)
+        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(<self|self>)"""
@@ -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
+        """Checks if the matrix product of itself with conjugate transpose is Identity UU+ = I"""
-          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]
         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],
+CNOT = TwoQubitOperator([[1, 0, 0, 0],
                          [0, 1, 0, 0],
                          [0, 0, 0, 1],
-                        [0, 0, 1, 0], ])
+                         [0, 0, 1, 0], ],
+                        TwoQubitPartial("C"),
+                        TwoQubitPartial("x"),
+                        I,
+                        X)
 
-
+C, x = CNOT.A, CNOT.B
-# 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")
 
 
 # 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):
+class QuantumCircuit(object):
-    HALT_STATE = "HALT"
-    RUNNING_STATE = "RUNNING"
-
     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)
+    qc = QuantumCircuit(3)
-    qp.add_row([H, C])
+    qc.add_row([H, C])
-    qp.add_row([_, _])
+    qc.add_row([_, _])
-    qp.add_row([_, x])
+    qc.add_row([_, x])
-    qp.print()
+    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)