diff --git a/lib_q_computer_math.py b/lib_q_computer_math.py index d6efe4f..2cc657b 100644 --- a/lib_q_computer_math.py +++ b/lib_q_computer_math.py @@ -1,5 +1,8 @@ from __future__ import annotations +import random +from collections import defaultdict +from functools import reduce from typing import Union import numpy as np @@ -34,6 +37,8 @@ class Matrix(object): def __eq__(self, other): if isinstance(other, Complex): return bool(np.allclose(self.m, other)) + elif isinstance(other, TwoQubitGate): + return False return np.allclose(self.m, other.m) def __mul_or_kron__(self, other): @@ -61,6 +66,14 @@ class Matrix(object): 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)) + + def x(self, other): + """Define outer product |0><0| looks like |0x0| which is 0.x(0)""" + return self.outer(other) + def _conjugate_transpose(self): return self.m.transpose().conjugate() @@ -68,10 +81,10 @@ class Matrix(object): return self.m.conjugate() def conjugate_transpose(self): - return State(self._conjugate_transpose()) + return self.__class__(self._conjugate_transpose()) def complex_conjugate(self): - return State(self._complex_conjugate()) + return self.__class__(self._complex_conjugate()) class State(Matrix): @@ -97,12 +110,18 @@ class State(Matrix): """If the inner (dot) product is zero, this vector is orthogonal to other""" return self | other == 0 - def measure(self, j): + def get_prob(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 + 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}" + choices = [format_str.format(i) for i in range(len(weights))] + return random.choices(choices, weights)[0] + class Operator(object): def __init__(self, func=None, *args, **kwargs): @@ -137,15 +156,88 @@ 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) + 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) + def __repr__(self): + if self.name: + return '-{}-'.format(self.name) + return str(self.m) + + +""" +Define States and Operators +""" + +_0 = State([[1], + [0]], + name='0') + +_1 = State([[0], + [1]], + name='1') + +_ = I = UnitaryOperator([[1, 0], + [0, 1]], + name="-") + +X = UnitaryOperator([[0, 1], + [1, 0]], + name="X") + +Y = UnitaryOperator([[0, -1j], + [1j, 0]], + name="Y") + +Z = UnitaryOperator([[1, 0], + [0, -1]], + name="Z") + +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) +# +# Decomposed CNOT : +# reverse engineered from +# https://quantumcomputing.stackexchange.com/questions/4252/how-to-derive-the-cnot-matrix-for-a-3-qbit-system-where-the-control-target-qbi +# +# CNOT(q1, I, q2): +# |0><0| x I_2 x I_2 + |1><1| x I_2 x X +# np.kron(np.kron(np.outer(_0.m, _0.m), np.eye(2)), np.eye(2)) + np.kron(np.kron(np.outer(_1.m, _1.m), np.eye(2)), X.m) +# +# +# CNOT(q1, q2): +# |0><0| x I + |1><1| x X +# 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 TwoQubitGate(object): + def __init__(self, rpr): + self.rpr = rpr + + def __repr__(self): + return str("-{}-".format(self.rpr)) + + +C = TwoQubitGate("C") +x = TwoQubitGate("x") + +# TODO: End Hacky way to define 2-qbit gate +########################################################### + 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 @@ -172,22 +264,16 @@ def test(): # 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]]) + # Pauli X gate flips the |0> to |1> and the |1> to |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], @@ -197,14 +283,14 @@ def test(): # 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 _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 - 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 _1.get_prob(0) == 0 # Probability for |1> in 0 is 0 + assert _1.get_prob(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 + 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 assert _0 * _0 == State([[1], [0], [0], [0]]) @@ -215,10 +301,6 @@ def test(): # 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 @@ -239,10 +321,21 @@ def test(): [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 + assert np.isclose(bell.get_prob(0b00), 0.5) # Probability for bell in 00 is 0.5 + assert np.isclose(bell.get_prob(0b01), 0) # Probability for bell in 01 is 0. + assert np.isclose(bell.get_prob(0b10), 0) # Probability for bell in 10 is 0. + 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... + assert _0.x(_0) == Matrix([[1, 0], + [0, 0]]) + assert _0.x(_1) == Matrix([[0, 1], + [0, 0]]) + assert _1.x(_0) == Matrix([[0, 0], + [1, 0]]) + assert _1.x(_1) == Matrix([[0, 0], + [0, 1]]) def naive_load_test(N): @@ -289,5 +382,135 @@ def naive_load_test(N): gc.collect() +class QuantumProcessor(object): + HALT_STATE = "HALT" + RUNNING_STATE = "RUNNING" + + def __init__(self, n_qubits: int): + self.n_qubits = n_qubits + self.steps = [[_0 for _ in range(n_qubits)], ] + self._called_add_row = 0 + self.current_step = 0 + self.current_quantum_state = None + self.current_state = self.HALT_STATE + self.reset() + + def add_row_step(self, row: int, step: int, qbit_state): + if len(self.steps) <= step: + self.steps += [[I for _ in range(self.n_qubits)] for _ in range(len(self.steps) - step + 1)] + self.steps[step][row] = qbit_state + + def add_step(self, step_data: list): + if len(step_data) != self.n_qubits: + raise RuntimeError("Length of step is: {}, should be: {}".format(len(step_data), self.n_qubits)) + step_i = len(step_data) + for row, qubit_state in enumerate(step_data): + self.add_row_step(row, step_i, qubit_state) + + def add_steps(self, steps_data: list): + for step_data in steps_data: + self.add_step(step_data) + + def add_row(self, row_data: list): + if self._called_add_row >= self.n_qubits: + raise RuntimeError("Adding more rows than qubits") + for step_i, qubit_state in enumerate(row_data): + self.add_row_step(self._called_add_row, step_i + 1, qubit_state) + self._called_add_row += 1 + + def add_rows(self, rows_data: list): + for row_data in rows_data: + 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) + return reduce((lambda x, y: x * y), step) + + def step(self): + if self.current_step == 0 and self.current_state == self.HALT_STATE: + self.current_state = self.RUNNING_STATE + if self.current_step >= len(self.steps): + self.current_state = self.HALT_STATE + raise RuntimeWarning("Halted") + running_step = self.steps[self.current_step] + step_quantum_state = self.compose_quantum_state(running_step) + 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_step += 1 + + def run(self): + for _ in self.steps: + self.step() + self.current_state = self.HALT_STATE + + def reset(self): + self.current_step = 0 + self.current_quantum_state = None + self.current_state = self.HALT_STATE + + def print(self): + print("=" * 3 * len(self.steps)) + for line_no in range(self.n_qubits): + line = '' + for step in self.steps: + state_repr = repr(step[line_no]) + line += state_repr + print(line) + print("=" * 3 * len(self.steps)) + + def measure(self): + if self.current_state != self.HALT_STATE: + raise RuntimeError("Processor is still running") + return self.current_quantum_state.measure() + + def run_n(self, n: int): + for i in range(n): + self.run() + result = self.measure() + print("Run {}: {}".format(i, result)) + self.reset() + + def get_sample(self, n: int): + rv = defaultdict(int) + for i in range(n): + self.run() + result = self.measure() + rv[result] += 1 + self.reset() + for k, v in sorted(rv.items(), key=lambda x: x[0]): + print("{}: {}".format(k, v)) + return rv + + +def test_quantum_processor(): + # Produce Bell state between 0 and 2 qubit + qp = QuantumProcessor(3) + qp.add_row([H, C]) + qp.add_row([_, _]) + qp.add_row([_, x]) + qp.print() + qp.get_sample(100) + + if __name__ == "__main__": - naive_load_test(23) + test() + test_quantum_processor()