auxilary methods for state

lib_q_computer_math.py

@@ -1,12 +1,11 @@
 import random
 from collections import defaultdict
 from functools import reduce
+from numbers import Complex
 from typing import Union
 
 import numpy as np
-from numbers import Complex
+import sympy as sp
-
-import sympy
 
 from load_test import sizeof_fmt
 
@@ -105,6 +104,9 @@ class Matrix(object):
     def complex_conjugate(self):
         return Matrix(self._complex_conjugate())
 
+    def __repr__(self):
+        return str(self.m)
+
 
 class Vector(Matrix):
     def __init__(self, m: ListOrNdarray = None, *args, **kwargs):
@@ -168,17 +170,29 @@ class State(Vector):
 
     def __repr__(self):
         if not self.name:
+            if np.count_nonzero(self.m == 1):
+                fmt_str = self.get_fmt_of_element()
+                index = np.where(self.m == [1])[0][0]
+                self.name = fmt_str.format(index)
+            else:
                 for well_known_state in well_known_states:
+                    try:
                         if self == well_known_state:
                             return repr(well_known_state)
+                    except:
+                        ...
                 for used_state in UNIVERSE_STATES:
+                    try:
                         if self == used_state:
                             return repr(used_state)
+                    except:
+                        ...
                 next_state_sub = ''.join([REPR_MATH_SUBSCRIPT_NUMBERS[int(d)] for d in str(len(UNIVERSE_STATES))])
                 self.name = '{}{}'.format(REPR_GREEK_PSI, next_state_sub)
         UNIVERSE_STATES.append(self)
-        matrix_rep = "{}".format(self.m).replace('[', '').replace(']', '').replace('\n', '|').strip()
+        # matrix_rep = "{}".format(self.m).replace('[', '').replace(']', '').replace('\n', '|').strip()
-        state_name = '|{}> = {}'.format(self.name, matrix_rep)
+        # state_name = '|{}> = {}'.format(self.name, matrix_rep)
+        state_name = '|{}>'.format(self.name)
         return state_name
 
     def norm(self):
@@ -199,12 +213,89 @@ class State(Vector):
         e_j = State([[1] if i == int(j) else [0] for i in range(len(self.m))])
         return np.absolute((e_j | self).m.item(0)) ** 2
 
+    def get_fmt_of_element(self):
+        return "{:0" + str(int(np.ceil(np.log2(len(self))))) + "b}"
+
     def measure(self):
         weights = [self.get_prob(j) for j in range(len(self))]
-        format_str = "{:0" + str(int(np.ceil(np.log2(len(weights))))) + "b}"
+        format_str = self.get_fmt_of_element()
         choices = [format_str.format(i) for i in range(len(weights))]
         return random.choices(choices, weights)[0]
 
+    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 pretty_print(self):
+        format_str = self.get_fmt_of_element() + " | {}"
+
+        for i, element in enumerate(self.m):
+            print(format_str.format(i, humanize_num(element[0])))
+
+    @classmethod
+    def from_string(cls, q):
+        try:
+            bin_repr = q[1:-1]
+            dec_bin = int(bin_repr, 2)
+            bin_length = 2 ** len(bin_repr)
+            arr = [[1] if i == dec_bin else [0] for i in range(bin_length)]
+        except:
+            raise Exception("State from string should be of the form |00..01> with all numbers either 0 or 1")
+        return cls(arr)
+
+    @staticmethod
+    def get_random_state_angles():
+        phi = np.random.uniform(0, 2 * np.pi)
+        t = np.random.uniform(0, 1)
+        theta = np.arccos(1 - 2 * t)
+        return theta, phi
+
+    def get_bloch_coordinates(self):
+        theta, phi = self.to_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):
+    """Helper method for creating state easily"""
+    if type(q) == str:
+        # e.g. |000>
+        if q[0] == '|' and q[-1] == '>':
+            return State.from_string(q)
+    elif type(q) == list:
+        # e.g. s([1,0]) => |0>
+        return State(np.reshape(q, (len(q), 1)))
+    return State(q)
+
+
+def humanize_num(fl, tolerance=1e-3):
+    if np.abs(fl) < tolerance:
+        return 0
+    if np.abs(fl.imag) < tolerance:
+        fl = fl.real
+    try:
+        return sp.nsimplify(fl, [sp.pi], tolerance, full=True)
+    except:
+        return sp.nsimplify(fl, [sp.pi], tolerance)
+
+
+def pp(m, tolerance=1e-3):
+    for element in m:
+        print(humanize_num(element, tolerance=tolerance))
+
+
+def humanize(m):
+    rv = []
+    for element in m:
+        rv.append(humanize(element))
+    return rv
+
 
 class Operator(object):
     def __init__(self, func=None, *args, **kwargs):
@@ -265,6 +356,18 @@ class UnitaryMatrix(SquareMatrix):
         return np.isclose(UU_, I).all()
 
 
+class HermitianMatrix(SquareMatrix):
+    def __init__(self, m: ListOrNdarray, *args, **kwargs):
+        """Represents a Hermitian matrix that satisfies U=U+"""
+        super().__init__(m, *args, **kwargs)
+        if not self._is_hermitian():
+            raise TypeError("Not a Hermitian matrix")
+
+    def _is_hermitian(self):
+        """Checks if its equal to the conjugate transpose U = U+"""
+        return np.isclose(self.m, self._conjugate_transpose()).all()
+
+
 class UnitaryOperator(LinearOperator, UnitaryMatrix):
     def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs):
         """UnitaryOperator inherits from both LinearOperator and a Unitary matrix
@@ -364,19 +467,7 @@ _m = State([[1 / np.sqrt(2)],
             [- 1 / np.sqrt(2)]],
            name='-')
 
-_00 = State([[1],
+well_known_states = [_p, _m]
-             [0],
-             [0],
-             [0]],
-            name='00')
-
-_11 = State([[0],
-             [0],
-             [0],
-             [1]],
-            name='11')
-
-well_known_states = [_0, _1, _p, _m]
 
 _ = I = UnitaryOperator([[1, 0],
                          [0, 1]],
@@ -414,6 +505,57 @@ C, x = CNOT.A, CNOT.B
 ###########################################################
 
 
+def assert_raises(exception, msg, callable, *args, **kwargs):
+    try:
+        callable(*args, **kwargs)
+    except exception as e:
+        assert str(e) == msg
+        return
+    assert False
+
+
+def assert_not_raises(exception, msg, callable, *args, **kwargs):
+    try:
+        callable(*args, **kwargs)
+    except exception as e:
+        if str(e) == msg:
+            assert False
+
+
+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?
+    h_not_u = [
+        [1, 0],
+        [0, 2],
+    ]
+    assert_not_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix, h_not_u)
+    assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, h_not_u)
+
+    u_not_h = [
+        [1, 0],
+        [0, 1j],
+    ]
+    assert_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix, u_not_h)
+    assert_not_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, u_not_h)
+
+    u_and_h = [
+        [0, 1],
+        [1, 0],
+    ]
+    assert_not_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix, u_and_h)
+    assert_not_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, u_and_h)
+
+    not_u_not_h = [
+        [1, 2],
+        [0, 1],
+    ]
+    assert_raises(TypeError, "Not a Hermitian matrix", HermitianMatrix, not_u_not_h)
+    assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, not_u_not_h)
+
+
 def test():
     # Test properties of Hilbert vector space
     # The four postulates of Quantum Mechanics
@@ -445,6 +587,9 @@ def test():
     assert _times_2.on(5) == 10
     assert _times_2(5) == 10
 
+    # Understanding the difference between unitary and hermitians
+    test_unitary_hermitian()
+
     # Pauli X gate flips the |0> to |1> and the |1> to |0>
     assert X | _1 == _0
     assert X | _0 == _1
@@ -456,8 +601,8 @@ def test():
                             [0]])
 
     # III: Measurement | A quantum measurement is described by an orthonormal basis |e_j>
-    #                    for state space. If the initial state of the system is |psi>
+    #                    for state space. If the initial state of the system is |ψ>
-    #                    then we get outcome j with probability pr(j) = |<e_j|psi>|^2
+    #                    then we get outcome j with probability pr(j) = |<e_j|ψ>|^2
 
     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
@@ -468,7 +613,11 @@ 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
 
-    # IV: Compositing | tensor/kronecker product when composing
+    # IV: Compositing | The state space of a composite physical system
+    #                   is the tensor product of the state spaces
+    #                   of the component physical systems.
+    #                   i.e. if we have systems numbered 1 through n, ψ_i,
+    #                   then the joint state is |ψ_1> ⊗ |ψ_2> ⊗ ... ⊗ |ψ_n>
     assert _0 * _0 == State([[1], [0], [0], [0]])
     assert _0 * _1 == State([[0], [1], [0], [0]])
     assert _1 * _0 == State([[0], [0], [1], [0]])
@@ -534,6 +683,16 @@ def test_to_from_angles():
         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
@@ -767,11 +926,6 @@ def test_quantum_processor():
 
 if __name__ == "__main__":
     test()
-    test_to_from_angles()
+    pp(_p.get_bloch_coordinates())
-    test_quantum_processor()
+    # test_to_from_angles()
-    # print(_1 | _0)
+    # test_quantum_processor()
-    # qubit = _00 + _11
-    # print(qubit)
-    # bell = CNOT.on(qubit)
-    # HI = I*H
-    # print(bell)