Merge remote-tracking branch 'origin/master'

06_krisis.py

@@ -29,7 +29,7 @@ def print_all_samples(name, all_samples):
     print("==============================================")
 
 
-def math_sim(q2func=State.from_angles, iterations=1000, sample_count=1):
+def math_sim(q2func=State.from_bloch_angles, iterations=1000, sample_count=1):
     all_samples = defaultdict(int)
     for i in range(iterations):
         # print("Running iteration {}".format(i))
@@ -45,7 +45,7 @@ def math_sim(q2func=State.from_angles, iterations=1000, sample_count=1):
         # print("theta: {:.2f} ({:.2f} deg) | phi: {:.2f} ({:.2f} deg)".format(
         #     theta, np.rad2deg(theta),
         #     phi, np.rad2deg(phi)))
-        q1 = State.from_angles(theta, phi)
+        q1 = State.from_bloch_angles(theta, phi)
         q2 = q2func(theta, phi)
 
         qc = QuantumCircuit(2, [[q1, q2], ])
@@ -131,7 +131,7 @@ def cirq_sim(q1func, q2func, iterations=1000, sample_count=1):
 
 
 def math_sim_0(*args, **kwargs):
-    math_sim(q2func=State.from_angles, *args, **kwargs)
+    math_sim(q2func=State.from_bloch_angles, *args, **kwargs)
 
 
 def math_sim_1(*args, **kwargs):

3sat.py

@@ -0,0 +1,106 @@
+import inspect
+import itertools
+import random
+from pprint import pprint
+
+from grover import classical_func_search_multi_rv
+
+
+class X(object):
+    """"""
+
+    def __init__(self, i):
+        self.i = i
+
+
+class SAT(object):
+    def __init__(self, predicates=3):
+        self.and_clauses = []
+        self.and_clauses_funcs = []
+        self.predicates = predicates
+
+    def or_(self, or_clauses):
+        def inner_or_(params):
+            for or_clause in or_clauses:
+                yes_no, var = or_clause[0], or_clause[1].i - 1
+                if yes_no(params[var]):
+                    return True
+            return False
+
+        return inner_or_
+
+    def and_(self, and_clauses):
+        def inner_and_(params):
+            for and_clause in and_clauses:
+                if not and_clause(params):
+                    return False
+            return True
+
+        return inner_and_
+
+    def yes(self, x):
+        return x
+
+    def no(self, x):
+        return not x
+
+    def generate(self, ands=1):
+        params = [X(i) for i in range(1, self.predicates + 1)]
+
+        def or_clause(or_clauses):
+            def inner_or_clause(params):
+                return self.or_(or_clauses)(params)
+
+            return inner_or_clause
+
+        for _ in range(ands):
+            or_clauses = list(zip(random.choices([self.yes, self.no], k=3), random.sample(params, k=3)))
+            self.and_clauses.append(or_clauses)
+            self.and_clauses_funcs.append(or_clause(or_clauses))
+
+        def inner_generate(args):
+            return self.and_(self.and_clauses_funcs)(args)
+
+        self.inner = inner_generate
+
+    def __call__(self, args):
+        return self.inner(args)
+
+    def show(self):
+        for j, _or_clause in enumerate(self.and_clauses):
+            for i, predicate in enumerate(_or_clause):
+                if i == 0:
+                    print("( ", end='')
+                trueness, var = predicate
+                if trueness == self.yes:
+                    print("x_{}".format(var.i), end='')
+                else:
+                    print("¬x_{}".format(var.i), end='')
+                if i + 1 != len(_or_clause):
+                    print(" ∨ ", end='')
+            if j + 1 != len(self.and_clauses):
+                print(" ) ∧ ")
+            else:
+                print(" )")
+
+
+def classical_3sat(func):
+    # Generate all possible true/false tupples for the 3-sat problem
+    input_range = list(itertools.product([True, False], repeat=func.predicates))
+    random.shuffle(input_range)
+    return classical_func_search_multi_rv(func, input_range)
+
+
+def main():
+    rv = []
+    for _ in range(100):
+        gen_3sat = SAT(predicates=4)
+        gen_3sat.generate(ands=8)
+        # gen_3sat.show()
+        sols = classical_3sat(gen_3sat)
+        rv.append(len(sols))
+    print(sorted(rv))
+
+
+if __name__ == "__main__":
+    main()

compiler.py

@@ -0,0 +1,9 @@
+from lib_q_computer_math import s
+
+
+def int_to_state(i):
+    return s("|{}>".format(bin(i)))
+
+
+if __name__ == "__main__":
+    print(int_to_state(42))

grover.py

@@ -0,0 +1,58 @@
+# http://twistedoakstudios.com/blog/Post2644_grovers-quantum-search-algorithm
+import itertools
+import random
+from pprint import pprint
+
+
+def classical_func_search_multi_rv(func, input_range):
+    """Grover’s algorithm takes a function, searches through
+    the implicit list of possible inputs to that function, and
+    returns inputs that causes the function to return true.
+    """
+    rv = []
+    for i, params in enumerate(input_range):
+        result = func(params)
+        if result:
+            rv.append(params)
+    return rv
+
+
+def classical_func_search_single_rv(func, input_range):
+    """Grover’s algorithm takes a function, searches through
+    the implicit list of possible inputs to that function, and
+    returns EXACTLY ONE input that causes the function to return true.
+    """
+    rv = None
+    for i, params in enumerate(input_range):
+        result = func(params)
+        if result and rv is None:
+            rv = result
+        else:
+            raise Exception("Exactly one result is needed")
+    return rv
+
+
+def _3sat(params):
+    # https://cstheory.stackexchange.com/questions/38538/oracle-construction-for-grovers-algorithm
+    # 3-SAT from here: https://qiskit.org/textbook/ch-applications/satisfiability-grover.html
+    x1, x2, x3 = params
+    return (not x1 or not x2 or not x3) and \
+           (x1 or not x2 or x3) and \
+           (x1 or x2 or not x3) and \
+           (x1 or not x2 or not x3) and \
+           (not x1 or x2 or x3)
+
+
+def classical_3sat(func):
+    # Generate all possible true/false tupples for the 3-sat problem
+    input_range = list(itertools.product([True, False], repeat=3))
+    random.shuffle(input_range)
+    return classical_func_search_multi_rv(func, input_range)
+
+
+def main():
+    pprint(classical_3sat(_3sat))
+
+
+if __name__ == "__main__":
+    main()

lib_q_computer_math.py

@@ -12,6 +12,7 @@ from load_test import sizeof_fmt
 ListOrNdarray = Union[list, np.ndarray]
 
 REPR_TENSOR_OP = "⊗"
+REPR_GREEK_BETA = "β"
 REPR_GREEK_PSI = "ψ"
 REPR_GREEK_PHI = "φ"
 REPR_MATH_KET = "⟩"
@@ -47,7 +48,7 @@ class Matrix(object):
     def __eq__(self, other):
         if isinstance(other, Complex):
             return bool(np.allclose(self.m, other))
-        elif isinstance(other, TwoQubitPartial):
+        elif isinstance(other, PartialQubit):
             return False
         return np.allclose(self.m, other.m)
 
@@ -104,10 +105,26 @@ class Matrix(object):
     def complex_conjugate(self):
         return Matrix(self._complex_conjugate())
 
+    @property
+    def human_m(self):
+        return humanize(self.m)
+
     def __repr__(self):
         return str(self.m)
 
 
+class HorizontalVector(Matrix):
+    """Horizontal vector is basically a list"""
+
+    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 len(self.m.shape) == 1
+
+
 class Vector(Matrix):
     def __init__(self, m: ListOrNdarray = None, *args, **kwargs):
         super().__init__(m, *args, **kwargs)
@@ -123,6 +140,7 @@ class State(Vector):
         """State vector representing quantum state"""
         super().__init__(m, *args, **kwargs)
         self.name = name
+        self.measurement_result = None
         if not self._is_normalized():
             raise TypeError("Not a normalized state vector")
 
@@ -138,16 +156,18 @@ class State(Vector):
         return theta, phi
 
     @classmethod
-    def from_angles(cls, theta, phi):
+    def from_bloch_angles(cls, theta, phi):
+        """Creates a state from angles in Bloch sphere"""
         theta, phi = cls._normalize_angles(theta, phi)
         m0 = np.cos(theta / 2)
         m1 = np.sin(theta / 2) * np.power(np.e, (1j * phi))
         m = m0 * _0 + m1 * _1
         return cls(m.m)
 
-    def to_angles(self):
+    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 be 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
@@ -168,6 +188,15 @@ class State(Vector):
         assert 0 <= phi <= 2 * np.pi
         return theta, phi
 
+    def rotate_x(self, theta):
+        return Rx(theta).on(self)
+
+    def rotate_y(self, theta):
+        return Ry(theta).on(self)
+
+    def rotate_z(self, theta):
+        return Rz(theta).on(self)
+
     def __repr__(self):
         if not self.name:
             if np.count_nonzero(self.m == 1):
@@ -217,18 +246,63 @@ class State(Vector):
         return "{:0" + str(int(np.ceil(np.log2(len(self))))) + "b}"
 
     def measure(self):
+        """
+        Measures in the computational basis.
+        Irreversable operation. Measuring again will result in the same result
+        TODO: Generalize the method so it takes a basis
+        TODO: Should we memoize per basis?
+              If it's measured twice, should it return the same state?
+              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?
+        :return: binary representation of the measured qubit (e.g. "011")
+        """
+        if self.measurement_result:
+            return self.measurement_result
         weights = [self.get_prob(j) for j in range(len(self))]
         format_str = self.get_fmt_of_element()
         choices = [format_str.format(i) for i in range(len(weights))]
-        return random.choices(choices, weights)[0]
+        self.measurement_result = random.choices(choices, weights)[0]
+        return self.measurement_result
 
-    def measure_n(self, n=1):
-        """measures n times"""
-        measurements = defaultdict(int)
-        for _ in range(n):
-            k = self.measure()
-            measurements[k] += 1
-        return measurements
+    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
+        """
+        m = mo.on(self) / np.sqrt(mo.get_prob(self))
+        return State(m)
+
+    def measure_partial(self, qubit_n):
+        """Partial measurement of state
+        Measures the n-th qubit with probability sum(a_n),
+        adjusting the probability of the rest of the state
+        """
+        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))
+        format_str = self.get_fmt_of_element()
+        # e.g. for state |000>:
+        # ['000', '001', '010', '011', '100', '101', '110', '111']
+        bin_repr = [format_str.format(i) for i in range(len(self))]
+        # e.g. for qbit_n = 2
+        # [0, 0, 1, 1, 0, 0, 1, 1]
+        partial_measurement_of_qbit = [int(b[qubit_n - 1]) for b in bin_repr]
+        # [0, 1, 4, 5]
+        indexes_for_p_0 = [i for i, index in enumerate(partial_measurement_of_qbit) if index == 0]
+        weights_0 = [self.get_prob(j) for j in indexes_for_p_0]
+        choices_0 = [format_str.format(i) for i in indexes_for_p_0]
+        measurement_result = random.choices(choices_0, weights_0)[0][qubit_n - 1]
+        # TODO: Verify if this is the correct collapse to lower dimension after partial measurement
+        #       https://www.youtube.com/watch?v=MG_9JWsrKtM&list=PL1826E60FD05B44E4&index=16
+        normalization_factor = np.sqrt(np.sum([self.m[i][0] for i in indexes_for_p_0]))
+        # TODO: This can be 0...
+        self.m = [
+            [(self.m[i][0] ** 2) / normalization_factor] for i in indexes_for_p_0
+        ]
+        return measurement_result
 
     def pretty_print(self):
         format_str = self.get_fmt_of_element() + " | {}"
@@ -255,30 +329,45 @@ class State(Vector):
         return theta, phi
 
     def get_bloch_coordinates(self):
-        theta, phi = self.to_angles()
+        theta, phi = self.to_bloch_angles()
         x = np.sin(theta) * np.cos(phi)
         y = np.sin(theta) * np.sin(phi)
         z = np.cos(theta)
         return [x, y, z]
 
 
-def s(q):
+def test_measure_partial():
+    state = s("|010>")
+    state.measure_partial(2)
+
+
+def normalize_state(state_vector: ListOrNdarray):
+    """Normalize a state by dividing by the square root of sum of the squares of states"""
+    norm_coef = np.sqrt(np.sum(np.array(state_vector) ** 2))
+    if norm_coef == 0:
+        raise TypeError("zero-sum vector")
+    return state_vector / norm_coef
+
+
+def s(q, name=None):
     """Helper method for creating state easily"""
     if type(q) == str:
         # e.g. |000>
         if q[0] == '|' and q[-1] == '>':
-            return State.from_string(q)
+            state = State.from_string(q)
+            state.name = name
+            return state
     elif type(q) == list:
         # e.g. s([1,0]) => |0>
-        return State(np.reshape(q, (len(q), 1)))
-    return State(q)
+        return State(np.reshape(q, (len(q), 1)), name=name)
+    return State(q, name=name)
 
 
 def humanize_num(fl, tolerance=1e-3):
-    if np.abs(fl) < tolerance:
-        return 0
     if np.abs(fl.imag) < tolerance:
         fl = fl.real
+    if np.abs(fl.real) < tolerance:
+        fl = 0 + 1j * fl.imag
     try:
         return sp.nsimplify(fl, [sp.pi], tolerance, full=True)
     except:
@@ -293,43 +382,56 @@ def pp(m, tolerance=1e-3):
 def humanize(m):
     rv = []
     for element in m:
-        rv.append(humanize(element))
+        if type(element) in [np.ndarray, list]:
+            rv.append(humanize(element))
+        else:
+            rv.append(humanize_num(element))
     return rv
 
 
 class Operator(object):
-    def __init__(self, func=None, *args, **kwargs):
+    def __init__(self, func: sp.Lambda, *args, **kwargs):
         """An Operator turns one function into another"""
         self.func = func
 
-    def on(self, *args, **kwargs):
-        return self(*args, **kwargs)
+    def on(self, *args):
+        return self(*args)
 
-    def __call__(self, *args, **kwargs):
-        return self.func(*args, **kwargs)
+    def __call__(self, *args):
+        return self.func(*args)
 
 
 class LinearOperator(Operator):
-    def __init__(self, func=None, *args, **kwargs):
+    def __init__(self, func: sp.Lambda, *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)
+        # 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
+        expr, vars_ = self.func.expr, self.func.variables
+        for x in vars_:
+            for y in vars_:
+                try:
+                    if not sp.Eq(sp.diff(expr, x, y), 0):
+                        return False
+                except TypeError:
+                    return False
         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
+
+
+def test_linear_operator():
+    # Operators turn one vector into another
+    # the times 2 operator should return the times two multiplication
+    _x = sp.Symbol('x')
+    _times_2 = LinearOperator(sp.Lambda(_x, 2 * _x))
+    assert _times_2.on(5) == 10
+    assert _times_2(5) == 10
+
+    assert_raises(TypeError, "Not a linear operator", LinearOperator, sp.Lambda(_x, _x ** 2))
 
 
 class SquareMatrix(Matrix):
@@ -342,6 +444,39 @@ class SquareMatrix(Matrix):
         return self.m.shape[0] == self.m.shape[1]
 
 
+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"""
+        self.name = name
+        Matrix.__init__(self, m=m, *args, **kwargs)
+        self.func = func or self.operator_func
+        LinearOperator.__init__(self, func=func, *args, **kwargs)
+
+    def _is_linear(self):
+        # 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):
+        return Vector(np.dot(self.m, other.m))
+
+    def get_eigens(self):
+        """ Returns (eigenvalue, eigenvector)
+        M|v> = λ|v>  ->
+          M   :is the operator (self)
+          |v> :is eigenstate of M
+          λ   :is the corresponding eigenvalue
+        """
+        eigenvalues, eigenvectors = np.linalg.eig(self.m)
+        rv = []
+        for i in range(0, len(eigenvectors)):
+            eigenvalue = eigenvalues[i]
+            eigenvector = HorizontalVector(eigenvectors[:, i])
+            rv.append((eigenvalue, eigenvector))
+        return rv
+
+
 class UnitaryMatrix(SquareMatrix):
     def __init__(self, m: ListOrNdarray, *args, **kwargs):
         """Represents a Unitary matrix that satisfies UU+ = I"""
@@ -368,13 +503,13 @@ class HermitianMatrix(SquareMatrix):
         return np.isclose(self.m, self._conjugate_transpose()).all()
 
 
-class UnitaryOperator(LinearOperator, UnitaryMatrix):
+class UnitaryOperator(LinearTransformation, UnitaryMatrix):
     def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs):
-        """UnitaryOperator inherits from both LinearOperator and a Unitary matrix
+        """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"""
-        self.name = name
         UnitaryMatrix.__init__(self, m=m, *args, **kwargs)
-        LinearOperator.__init__(self, func=self.operator_func, *args, **kwargs)
+        LinearTransformation.__init__(self, m=m, func=self.operator_func, *args, **kwargs)
+        self.name = name
 
     def operator_func(self, other):
         return State(np.dot(self.m, other.m))
@@ -385,8 +520,26 @@ class UnitaryOperator(LinearOperator, 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)
+class HermitianOperator(LinearTransformation, HermitianMatrix):
+    def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs):
+        """HermitianMatrix 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"""
+        return np.dot(self.m, other.m)
+
+    def __repr__(self):
+        if self.name:
+            return '-{}-'.format(self.name)
+        return str(self.m)
+
+
+# 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
@@ -402,7 +555,7 @@ class UnitaryOperator(LinearOperator, UnitaryMatrix):
 # 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):
+class PartialQubit(object):
     def __init__(self, rpr):
         self.rpr = rpr
         self.operator = None
@@ -411,12 +564,12 @@ class TwoQubitPartial(object):
         return str("-{}-".format(self.rpr))
 
 
-C_ = TwoQubitPartial("C")
-x_ = TwoQubitPartial("x")
+C_partial = PartialQubit("C")
+x_partial = PartialQubit("x")
 
 
 class TwoQubitOperator(UnitaryOperator):
-    def __init__(self, m: ListOrNdarray, A: TwoQubitPartial, B: TwoQubitPartial,
+    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
@@ -425,13 +578,16 @@ class TwoQubitOperator(UnitaryOperator):
         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 CONTROL and TARGET need to be defined in the same step exactly once")
+            raise RuntimeError("Both {} and {} need to be defined in the same step exactly once".format(
+                self.A, self.B
+            ))
 
     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:
@@ -459,6 +615,7 @@ _1 = State([[0],
             [1]],
            name='1')
 
+# |+> and |-> states
 _p = State([[1 / np.sqrt(2)],
             [1 / np.sqrt(2)]],
            name='+')
@@ -467,7 +624,32 @@ _m = State([[1 / np.sqrt(2)],
             [- 1 / np.sqrt(2)]],
            name='-')
 
-well_known_states = [_p, _m]
+# 4 Bell states
+b_00 = State([[1 / np.sqrt(2)],
+              [0],
+              [0],
+              [1 / np.sqrt(2)]],
+             name='B00')
+
+b_01 = State([[0],
+              [1 / np.sqrt(2)],
+              [1 / np.sqrt(2)],
+              [0]],
+             name='B01')
+
+b_10 = State([[1 / np.sqrt(2)],
+              [0],
+              [0],
+              [-1 / np.sqrt(2)]],
+             name='B10')
+
+b_11 = State([[0],
+              [1 / np.sqrt(2)],
+              [-1 / np.sqrt(2)],
+              [0]],
+             name='B11')
+
+well_known_states = [_p, _m, b_00, b_01, b_10, b_11]
 
 _ = I = UnitaryOperator([[1, 0],
                          [0, 1]],
@@ -485,24 +667,53 @@ Z = UnitaryOperator([[1, 0],
                      [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
+#       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
+#
+
+def Rx(theta):
+    return UnitaryOperator([[np.cos(theta / 2), -1j * np.sin(theta / 2)],
+                            [-1j * np.sin(theta / 2), np.cos(theta / 2)]],
+                           name="Rx")
+
+
+def Ry(theta):
+    return UnitaryOperator([[np.cos(theta / 2), -np.sin(theta / 2)],
+                            [np.sin(theta / 2), np.cos(theta / 2)]],
+                           name="Ry")
+
+
+def Rz(theta):
+    return UnitaryOperator([[np.power(np.e, -1j * theta / 2), 0],
+                            [0, np.power(np.e, 1j * theta / 2)]],
+                           name="Rz")
+
+
 H = UnitaryOperator([[1 / np.sqrt(2), 1 / np.sqrt(2)],
                      [1 / np.sqrt(2), -1 / np.sqrt(2)], ],
                     name="H")
 
-CNOT = TwoQubitOperator([[1, 0, 0, 0],
-                         [0, 1, 0, 0],
-                         [0, 0, 0, 1],
-                         [0, 0, 1, 0], ],
-                        TwoQubitPartial("C"),
-                        TwoQubitPartial("x"),
-                        I,
-                        X)
+CNOT = TwoQubitOperator([
+    [1, 0, 0, 0],
+    [0, 1, 0, 0],
+    [0, 0, 0, 1],
+    [0, 0, 1, 0],
+], C_partial, x_partial, I, X)
 
-C, x = CNOT.A, CNOT.B
-
-
-# TODO: End Hacky way to define 2-qbit gate
-###########################################################
+# 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)
 
 
 def assert_raises(exception, msg, callable, *args, **kwargs):
@@ -525,8 +736,8 @@ def assert_not_raises(exception, msg, callable, *args, **kwargs):
 def test_unitary_hermitian():
     # Unitary is UU+ = I; Hermitian is U = U+
     # Matrixes could be either, neither or both
-    # Quantum operators are described by unitary transformations
-    # TODO: What are Hermitians?
+    # 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
     h_not_u = [
         [1, 0],
         [0, 2],
@@ -556,6 +767,64 @@ def test_unitary_hermitian():
     assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, not_u_not_h)
 
 
+class MeasurementOperator(HermitianOperator):
+    """Measurement operators are Hermitians: <ψ|M†_m M_m|ψ>"""
+
+    @classmethod
+    def create_from_prob(cls, matrix: Matrix, *args, **kwargs):
+        """returns |M†_m><M_m|"""
+        return cls(matrix.conjugate_transpose().x(matrix).m, *args, **kwargs)
+
+    def get_prob(self, state: State):
+        """Returns result of <ψ|M†_m M_m|ψ>
+        """
+        # <ψ|
+        state_ct = state.conjugate_transpose()
+        #  This is: <ψ|   .   M†_m M_m . |ψ>
+        return state_ct.m.dot(self.m.dot(state.m)).item()
+
+
+def test_measurement_ops():
+    m0 = MeasurementOperator.create_from_prob(Matrix([1, 0]))
+    m1 = MeasurementOperator.create_from_prob(Matrix([0, 1]))
+    assert m0 == Matrix([[1, 0],
+                         [0, 0]])
+    assert m1 == Matrix([[0, 0],
+                         [0, 1]])
+
+    # p(0) -> probability of measurement to yield a 0
+    assert m0.get_prob(_0) == 1.0
+    assert m1.get_prob(_0) == 0.0
+    assert m0.get_prob(_1) == 0.0
+    assert m1.get_prob(_1) == 1.0
+
+    # Post-state measurement of qubit with operator
+    assert _p.measure_with_op(m0) == _0
+    assert _p.measure_with_op(m1) == _1
+    assert _m.measure_with_op(m0) == _0
+    assert _m.measure_with_op(m1) == s([0, -1])
+
+
+def abs_squared(x):
+    return np.abs(x) ** 2
+
+
+def testRotPauli():
+    # From http://www.vcpc.univie.ac.at/~ian/hotlist/qc/talks/bloch-sphere.pdf
+    # / slide 11:
+    # However, the only measurable quantities are the probabilities
+    # |α|**2 and |β|**2, so multiplying the state by an arbitrary factor
+    # exp(iγ) (a global phase) has no observable consequences
+    assert np.allclose(abs_squared(Rx(np.pi).m), abs_squared(X.m))
+    assert np.allclose(abs_squared(Ry(np.pi).m), abs_squared(Y.m))
+    assert np.allclose(abs_squared(Rz(np.pi).m), abs_squared(Z.m))
+
+
+def test_eigenstuff():
+    assert LinearTransformation(m=[[1, 0], [0, 0]]).get_eigens() == \
+           [(1.0, HorizontalVector([1., 0.])), (0., HorizontalVector([0., 1.]))]
+
+
 def test():
     # Test properties of Hilbert vector space
     # The four postulates of Quantum Mechanics
@@ -579,13 +848,13 @@ def test():
 
     assert _0 | _1 == (_1 | _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
     #
-    # Operators turn one vector into another
-    # the times 2 operator should return the times two multiplication
-    _times_2 = Operator(lambda x: 2 * x)
-    assert _times_2.on(5) == 10
-    assert _times_2(5) == 10
+    test_linear_operator()
+
+    test_eigenstuff()
 
     # Understanding the difference between unitary and hermitians
     test_unitary_hermitian()
@@ -600,9 +869,18 @@ def test():
     assert Y | _1 == State([[-1j],
                             [0]])
 
+    # 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) = |<e_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
+    #       are immeasurable using 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
     assert _0.get_prob(1) == 0  # Probability for |0> in 1 is 0
@@ -613,6 +891,9 @@ def test():
     assert np.isclose(_p.get_prob(0), 0.5)  # Probability for |+> in 0 is 0.5
     assert np.isclose(_p.get_prob(1), 0.5)  # Probability for |+> in 1 is 0.5
 
+    # Test measurement operators
+    test_measurement_ops()
+
     # IV: Compositing | The state space of a composite physical system
     #                   is the tensor product of the state spaces
     #                   of the component physical systems.
@@ -652,7 +933,7 @@ def test():
     assert np.isclose(bell.get_prob(0b11), 0.5)  # Probability for bell in 11 is 0.5
 
     ################################
-    # TODO: Don't know where outer product fits...
+    # TODO: Don't know where outer product fits - something about density operator?
     assert _0.x(_0) == Matrix([[1, 0],
                                [0, 0]])
     assert _0.x(_1) == Matrix([[0, 1],
@@ -665,34 +946,24 @@ def test():
 
 def test_to_from_angles():
     for q in [_0, _1, _p, _m]:
-        angles = q.to_angles()
-        s = State.from_angles(*angles)
+        angles = q.to_bloch_angles()
+        s = State.from_bloch_angles(*angles)
         assert q == s
 
-    assert State.from_angles(0, 0) == _0
-    assert State.from_angles(np.pi, 0) == _1
-    assert State.from_angles(np.pi / 2, 0) == _p
-    assert State.from_angles(np.pi / 2, np.pi) == _m
+    assert State.from_bloch_angles(0, 0) == _0
+    assert State.from_bloch_angles(np.pi, 0) == _1
+    assert State.from_bloch_angles(np.pi / 2, 0) == _p
+    assert State.from_bloch_angles(np.pi / 2, np.pi) == _m
     for theta, phi, qbit in [
         (0, 0, _0),
         (np.pi, 0, _1),
         (np.pi / 2, 0, _p),
         (np.pi / 2, np.pi, _m),
     ]:
-        s = State.from_angles(theta, phi)
+        s = State.from_bloch_angles(theta, phi)
         assert s == qbit
 
 
-def test_measure_n():
-    qq = State([
-        [0.5],
-        [0.5],
-        [0.5],
-        [0.5],
-    ])
-    qq.measure_n(100)
-
-
 def naive_load_test(N):
     import os
     import psutil
@@ -778,24 +1049,12 @@ class QuantumCircuit(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)
+        partials = [op for op in step if isinstance(op, PartialQubit)]
+        # 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()
+            return two_qubit_op.compose(step)
         return reduce((lambda x, y: x * y), step)
 
     def step(self):
@@ -849,8 +1108,8 @@ class QuantumProcessor(object):
         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)
+        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)
@@ -916,16 +1175,38 @@ class QuantumProcessor(object):
 def test_quantum_processor():
     # Produce Bell state between 0 and 2 qubit
     qc = QuantumCircuit(3)
-    qc.add_row([H, C])
+    qc.add_row([H, C_partial])
     qc.add_row([_, _])
-    qc.add_row([_, x])
+    qc.add_row([_, x_partial])
     qc.print()
     qp = QuantumProcessor(qc)
-    qp.get_sample(100)
+    qp.print_sample(qp.get_sample(100))
+
+
+def test_light():
+    # TODO: Are these measurement operators the correct way to represent hor/ver/diag filter?
+    #       No, becuase they are not unitaries
+    hor_filter = MeasurementOperator.create_from_prob(Matrix(_0.m), name='h')
+    diag_filter = MeasurementOperator.create_from_prob(Matrix(_p.m), name='d')
+    ver_filter = MeasurementOperator.create_from_prob(Matrix(_1.m), name='v')
+
+    # create random light/photon
+    random_pol = [[np.random.uniform(0, 1)], [np.random.uniform(0, 1)]]
+    random_light = State(normalize_state(random_pol))
+    qc = QuantumCircuit(1, initial_steps=[[random_light]])
+
+    # add three filters as operators
+    # qc.add_row([hor_filter, ver_filter])
+    qc.add_row([hor_filter, diag_filter, ver_filter])
+    qc.print()
+    qp = QuantumProcessor(qc)
+    # TODO: This doesn't work - dealing with non-normalized state vectors
+    #       - When measured in horizontal filter, the state vector has a positive |0> and 0 for |1>
+    qp.print_sample(qp.get_sample(100))
 
 
 if __name__ == "__main__":
-    test()
-    pp(_p.get_bloch_coordinates())
-    # test_to_from_angles()
+    # test()
     # test_quantum_processor()
+    # test_light()
+    test_measure_partial()

requirements.txt

@@ -9,7 +9,7 @@ cryptography==2.8
 cvxopt==1.2.4
 cycler==0.10.0
 Cython==0.29.14
-dataclasses==0.7
+dataclasses==0.6
 decorator==4.4.0
 dill==0.3.1.1
 dlx==1.0.4