compose two qubit gate

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: <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 +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
+        """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):
+        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")
+C_ = TwoQubitGate("C")
-x = TwoQubitGate("x")
+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 any([type(s) is TwoQubitGate for s in step]):
-            if not (step.count(C) == 1 and step.count(x) == 1):
+            two_qubit_gates = filter(lambda s: type(s) is TwoQubitGate, step)
-                raise RuntimeError("Both CONTROL and TARGET need to be defined in the same step exactly once")
+            state = []
-            # TODO: Hacky way to do CNOT
+            for two_qubit_gate in two_qubit_gates:
-            #       Should generalize for a 2-Qubit gate
+                two_qubit_gate.operator.verify_step(step)
-            # _0.x(_0) * Matrix(I.m) + _1.x(_1) * Matrix(X.m)
+                state = two_qubit_gate.operator.compose(step, state)
-            outer_0, outer_1 = [], []
+            return state
-            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)
         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):