from __future__ import annotations from typing import Union import numpy as np from numbers import Complex from load_test import sizeof_fmt ListOrNdarray = Union[list, np.ndarray] class Matrix(object): """Represents a quantum state as a vector in Hilbert space""" def __init__(self, m: ListOrNdarray = None, *args, **kwargs): """ Can be initialized with a matrix, e.g. for |0> this is [[0],[1]] :param m: a matrix representing the quantum state :param name: the name of the state (optional) """ if m is None: self.m = np.array([]) elif type(m) is list: self.m = np.array(m) elif type(m) is np.ndarray: self.m = m else: raise TypeError("m needs to be a list or ndarray type") def __add__(self, other): return other.__class__(self.m + other.m) def __eq__(self, other): if isinstance(other, Complex): return bool(np.allclose(self.m, other)) return np.allclose(self.m, other.m) def __mul_or_kron__(self, other): """ Multiplication with a number is Linear op; with another state is a composition via the Kronecker/Tensor product """ if isinstance(other, Complex): return self.__class__(self.m * other) elif isinstance(other, Matrix): return other.__class__(np.kron(self.m, other.m)) raise NotImplementedError def __rmul__(self, other): return self.__mul_or_kron__(other) def __mul__(self, other): return self.__mul_or_kron__(other) def __or__(self, other): """Define inner product: """ return other.__class__(np.dot(self._conjugate_transpose(), other.m)) def __len__(self): return len(self.m) def _conjugate_transpose(self): return self.m.transpose().conjugate() def _complex_conjugate(self): return self.m.conjugate() def conjugate_transpose(self): return State(self._conjugate_transpose()) def complex_conjugate(self): return State(self._complex_conjugate()) class State(Matrix): def __init__(self, m: ListOrNdarray = None, name: str = '', *args, **kwargs): """An Operator turns one function into another""" super().__init__(m) self.name = name def __repr__(self): if self.name: return '|{}>'.format(self.name) return str(self.m) def norm(self): """Norm/Length of the vector = sqrt()""" return self.length() def length(self): """Norm/Length of the vector = sqrt()""" return np.sqrt((self | self).m).item(0) def is_orthogonal(self, other): """If the inner (dot) product is zero, this vector is orthogonal to other""" return self | other == 0 def measure(self, j): """pr(j) = ||^2""" # j_th basis 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 class Operator(object): def __init__(self, func=None, *args, **kwargs): """An Operator turns one function into another""" self.func = func def on(self, *args, **kwargs): return self(*args, **kwargs) def __call__(self, *args, **kwargs): return self.func(*args, **kwargs) class UnitaryMatrix(Matrix): def __init__(self, m: ListOrNdarray, *args, **kwargs): """Represents a Unitary matrix that satisfies UU+ = I""" super().__init__(m, *args, **kwargs) if not self._is_unitary(): 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] UU_ = np.dot(self._conjugate_transpose(), self.m) I = np.eye(self.m.shape[0]) return is_square and np.isclose(UU_, I).all() class UnitaryOperator(Operator, UnitaryMatrix): def __init__(self, m: ListOrNdarray, name: str = '', *args, **kwargs): """UnitaryOperator inherits from both Operator and a Unitary matrix It is used to act on a State vector by defining the operator to be the dot product""" UnitaryMatrix.__init__(self, m=m, name=name, *args, **kwargs) Operator.__init__(self, func=lambda other: State(np.dot(self.m, other.m)), *args, **kwargs) def test(): _0 = State([[1], [0]], name='0') _1 = State([[0], [1]], name='1') _p = State([[1 / np.sqrt(2)], [1 / np.sqrt(2)]], name='+') m = 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 # known as the state space of the system. If the system is closed # then the system is described completely by its state vector # which is a unit vector in the space. # Mathematically, this vector space is also a function space assert _0 + _1 == _1 + _0 # commutativity of vector addition assert _0 + (_1 + _p) == (_0 + _1) + _p # associativity of vector addition assert 8 * (_0 + _1) == 8 * _0 + 8 * _1 # Linear when multiplying by constants assert _0 | _0 == 1 # parallel have 1 product assert _0 | _1 == 0 # orthogonal have 0 product assert _0.is_orthogonal(_1) assert _1 | (8 * _0) == 8 * (_1 | _0) # Inner product is linear multiplied by constants assert _p | (_1 + _0) == (_p | _1) + (_p | _0) # Inner product is linear in superpos of vectors assert np.isclose(_1.length(), 1.0) # all of the vector lengths are normalized assert np.isclose(_0.length(), 1.0) assert np.isclose(_p.length(), 1.0) assert _0 | _1 == (_1 | _0).complex_conjugate() # non-commutative inner product # II: Dynamics | The evolution of a closed system is described by a unitary transformation # # Operators turn one vector into another # Pauli X gate flips the |0> to |1> and the |1> to |0> # 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 X = UnitaryOperator([[0, 1], [1, 0]]) assert X | _1 == _0 assert X | _0 == _1 # Test the Y Pauli operator with complex number literal notation Y = UnitaryOperator([[0, -1j], [1j, 0], ]) assert Y | _0 == State([[0], [1j]]) assert Y | _1 == State([[-1j], [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> # then we get outcome j with probability pr(j) = ||^2 assert _0.measure(0) == 1 # Probability for |0> in 0 is 1 assert _0.measure(1) == 0 # Probability for |0> in 1 is 0 assert _1.measure(0) == 0 # Probability for |1> in 0 is 0 assert _1.measure(1) == 1 # Probability for |1> in 1 is 1 assert np.isclose(_p.measure(0), 0.5) # Probability for |+> in 0 is 0.5 assert np.isclose(_p.measure(1), 0.5) # Probability for |+> in 1 is 0.5 # IV: Compositing | tensor/kronecker product when composing assert _0 * _0 == State([[1], [0], [0], [0]]) assert _0 * _1 == State([[0], [1], [0], [0]]) assert _1 * _0 == State([[0], [0], [1], [0]]) assert _1 * _1 == State([[0], [0], [0], [1]]) # CNOT applies a control qubit on a target. # If the control is a |0>, target remains unchanged. # If the control is a |1>, target is flipped. CNOT = UnitaryOperator([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], ]) assert CNOT.on(_0 * _0) == _0 * _0 assert CNOT.on(_0 * _1) == _0 * _1 assert CNOT.on(_1 * _0) == _1 * _1 assert CNOT.on(_1 * _1) == _1 * _0 # ALL FOUR NOW - Create a Bell state (I) that has H|0> (II), measures (III) and composition in CNOT (IV) # Bell state # First - create a superposition H = UnitaryOperator([[1 / np.sqrt(2), 1 / np.sqrt(2)], [1 / np.sqrt(2), -1 / np.sqrt(2)], ]) superpos = H | _0 assert superpos == _p # Then CNOT the superposition with a |0> qubit bell = CNOT | (superpos * _0) assert bell == State([[1 / np.sqrt(2)], [0.], [0.], [1 / np.sqrt(2)], ]) assert np.isclose(bell.measure(0b00), 0.5) # Probability for bell in 00 is 0.5 assert np.isclose(bell.measure(0b01), 0) # Probability for bell in 01 is 0. assert np.isclose(bell.measure(0b10), 0) # Probability for bell in 10 is 0. assert np.isclose(bell.measure(0b11), 0.5) # Probability for bell in 11 is 0.5 def naive_load_test(N): import os import psutil import gc from time import time from sys import getsizeof print("{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format( "qbits", "kron_len", "mem_used", "mem_per_q", "getsizeof", "getsiz/len", "nbytes", "nbytes/len", "time")) _0 = State([[1], [0]], name='0') process = psutil.Process(os.getpid()) mem_init = process.memory_info().rss for i in range(2, N + 1): start = time() m = _0 for _ in range(i): m = m * _0 len_m = len(m) elapsed = time() - start mem_b = process.memory_info().rss - mem_init print("{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format( i, len(m), sizeof_fmt(mem_b), sizeof_fmt(mem_b / len_m), sizeof_fmt(getsizeof(m)), sizeof_fmt(getsizeof(m) / len_m), sizeof_fmt(m.m.nbytes), sizeof_fmt(m.m.nbytes / len_m), np.round(elapsed, 2))) gc.collect() if __name__ == "__main__": naive_load_test(23)