diff --git a/06_krisis.py b/06_krisis.py index 75dc7e6..7ffa6a9 100644 --- a/06_krisis.py +++ b/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): diff --git a/3sat.py b/3sat.py new file mode 100644 index 0000000..d39d56f --- /dev/null +++ b/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() diff --git a/compiler.py b/compiler.py new file mode 100644 index 0000000..4fc9e2a --- /dev/null +++ b/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)) diff --git a/grover.py b/grover.py new file mode 100644 index 0000000..986ffe4 --- /dev/null +++ b/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() diff --git a/lib_q_computer_math.py b/lib_q_computer_math.py index bb7eca4..2ea8b79 100644 --- a/lib_q_computer_math.py +++ b/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> + """ + # <ψ| + 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) = ||^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() diff --git a/papers/bloch-sphere.pdf b/papers/bloch-sphere.pdf new file mode 100644 index 0000000..aa5e710 Binary files /dev/null and b/papers/bloch-sphere.pdf differ diff --git a/requirements.txt b/requirements.txt index c014b97..366f029 100644 --- a/requirements.txt +++ b/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