Unverified Commit ad4637d6 authored by dwuggh's avatar dwuggh
Browse files

import

parents
__pycache__/
\ No newline at end of file
import numpy as np
from functools import reduce
from utils import *
from QOperator import *
from QChannel import *
class DensityOperator(QOperator):
def __init__(self, qubits, mat = None):
qnum = np.max(qubits)
dim = 2 ** qnum
# qubits = [i for i in range(qnum)]
# ground state: | 0000...0 >
# generate ground state density matrix
if mat is None:
mat = np.zeros((dim, dim))
mat[0, 0] = 1
# the following should be equivalent:
# ground_state = np.zeros(dim)
# ground_state[0] = 1
# mat = np.outer(ground_state, ground_state)
super().__init__(qubits, mat)
def print(self):
print("qubits: ", self.qubits)
print("shape: ", self.operator.shape)
print(self.operator)
def merge(self, other):
op = multiply(self, other)
self.operator = op.operator
self.qubits = op.qubits
def evolution(self, unitary: QOperator, p_g = 0, on_qubits = None):
# ρ = U ρ U^\dagger
if on_qubits is not None:
unitary.alter_qubits(on_qubits)
result = multiply(unitary, multiply(self, unitary.dagger()))
# big_unitary = unitary.broadcast(self.qnum())
# self.operator = np.matmul(np.matmul(big_unitary, self.operator), big_unitary.transpose())
self.operator = result.operator
self.qubits = result.qubits
if p_g != 0:
self.channel(depolarizing_channel(p_g, unitary.qubits))
return result
# TODO rename this function
def channel(self, channel: QChannel):
# qmax = np.max([self.qmax(), channel.qmax])
# result = sum(map(lambda E: np.matmul(np.matmul(E.operator, self.operator), E.dagger()), channel.kraus_operators))
result = reduce(add, (multiply(multiply(E, self), E.dagger()) for E in channel.kraus_operators))
self.operator = result.operator
self.qubits = result.qubits
return result
# bell measurement: only result 00 and 11 will be reserved
def bell_meausure(self, q1, q2, pauli, p_m = 0):
channel1 = measure_x(q1) if pauli == 'x' else measure_z(q2)
channel2 = measure_x(q1) if pauli == 'x' else measure_z(q2)
P_00 = multiply(channel1.kraus_operators[0], channel2.kraus_operators[0])
P_11 = multiply(channel1.kraus_operators[1], channel2.kraus_operators[1])
# this channel will cause probability loss because of post-selection
loss_channel = QChannel([P_00, P_11])
self.print()
self.channel(loss_channel)
self.print()
operator = self.partial_trace([q1, q2])
self.operator = operator.operator
self.qubits = operator.qubits
self.print()
# measurement channel
# E_0: measurement result 1, E_1: measurement result -1
def measure_x(q = 0):
E_0 = 0.5 * np.array([[1, 1], [1, 1]], dtype=np.complex64)
E_1 = 0.5 * np.array([[1, -1], [-1, 1]], dtype=np.complex64)
return QChannel([QOperator([q], E_0), QOperator([q], E_1)])
def measure_z(q = 0):
E_0 = np.array([[1, 0], [0, 0]], dtype=np.complex64)
E_1 = np.array([[0, 0], [0, 1]], dtype=np.complex64)
return QChannel([QOperator([q], E_0), QOperator([q], E_1)])
# generate bell pair in werner form
def bell_pair(p_n, qubits = [0, 1]):
Φ = np.array([1, 0, 0, 1], dtype=np.complex64) / 2 ** 0.5
mat = (1 - 4 / 3 * p_n) * np.outer(Φ, Φ) + (p_n / 3) * np.identity(4)
return DensityOperator(qubits, mat)
import numpy as np
from functools import reduce
from utils import *
from QOperator import *
class QChannel(object):
def __init__(self, kraus_operators):
# array of QOperators
self.kraus_operators = kraus_operators
# collection of all qubits
self.qubits = []
for op in self.kraus_operators:
for q in op.qubits:
if not contains(self.qubits, q):
self.qubits.append(q)
self.qubits = np.array(self.qubits, dtype=np.int32)
self.qmax = np.max(self.qubits) + 1
# broadcast every operator
# for op in self.kraus_operators:
# op = op.broadcast(qnum)
def depolarizing_channel(p, qubits = [0]) -> QChannel:
operators = []
qnum = len(qubits)
dim = 2 ** qnum
# (1 - p) I ρ I: E_0 = \sqrt{1 - p} I
E_0 = np.sqrt(1 - p) * np.identity(dim)
operators.append(QOperator(qubits, E_0))
# 4 ** qnum - 1 other operators
for i in range(4 ** qnum):
# i = 0 is just the identity operator, which has already been considered above
if i == 0:
continue
indices = np.flip(get_n_digits(i, 4, qnum))
# qi: a tuple, (q, i)
ps = map(lambda qi: pauli(qi[1], qi[0]), list(enumerate(indices)))
op = reduce(multiply, ps)
op.operator *= np.sqrt(p / (4 ** qnum - 1))
operators.append(op)
return QChannel(operators)
import numpy as np
from utils import *
class QOperator(object):
def __init__(self, qubits, operator: np.ndarray):
self.qubits = np.array(qubits)
self.operator = np.array(operator, dtype=np.complex64)
def qnum(self):
return self.qubits.size
def qmax(self):
return np.max(self.qubits) + 1
def dim(self):
# the following should be equivalent:
# return 2 ** self.qubits.size
return self.operator.shape[0]
def alter_qubits(self, qubits):
self.qubits = np.array(qubits)
# construct gate in the big hilbert space
# G = I ⊗ G
def broadcast(self, qnum: int):
dim = 2 ** qnum
big_operator = np.zeros((dim, dim), dtype=np.complex64)
for i in range(dim):
for j in range(dim):
di = np.flip(get_bin_digits(i, qnum))
dj = np.flip(get_bin_digits(j, qnum))
# coordinates in the small HIlbert space
gi = from_bin_digits(np.flip(di[self.qubits]))
gj = from_bin_digits(np.flip(dj[self.qubits]))
# check for identity
id = 1
for k in range(len(di)):
if id == 1 and not contains(self.qubits, k) and di[k] != dj[k]:
id = 0
val = np.complex64(id) * self.operator[gi, gj]
big_operator[i, j] = val
return big_operator
# def rearrange(self, qubits_new):
# pass
# get adjoint operator
def dagger(self):
t = self.operator.transpose()
t = t.conj()
return QOperator(self.qubits, t)
def partial_trace(self, traced_qubits):
dim = self.dim()
qnum = self.qnum()
qubits_reserved = np.fromiter((x for x in self.qubits if not contains(traced_qubits, x)), np.int32)
dim_new = dim - 2 ** len(traced_qubits)
# for simplicity, we treat scalar as a matrix of dim (1, 1)
dim_new = dim_new if dim_new > 0 else 1
# the traced matrix
mat_new = np.zeros((dim_new, dim_new), dtype=np.complex64)
for i in range(dim):
for j in range(dim):
di = np.flip(get_bin_digits(i, qnum))
dj = np.flip(get_bin_digits(j, qnum))
tr_i = from_bin_digits(np.flip(di[traced_qubits]))
tr_j = from_bin_digits(np.flip(dj[traced_qubits]))
i_new = from_bin_digits(np.flip(di[qubits_reserved]))
j_new = from_bin_digits(np.flip(dj[qubits_reserved]))
if tr_i == tr_j:
# print(i, j, i_new, j_new)
mat_new[i_new, j_new] += self.operator[i, j]
return QOperator(qubits_reserved, mat_new)
# multiply 2 QOperators, often needs to broadcast them
def multiply(lhs: QOperator, rhs: QOperator):
q1 = lhs.qubits
q2 = rhs.qubits
qmax = np.max(np.concatenate((q1, q2))) + 1
lhs_big = lhs.broadcast(qmax)
rhs_big = rhs.broadcast(qmax)
result_mat = np.matmul(lhs_big, rhs_big)
result_qubits = [i for i in range(qmax)]
return QOperator(result_qubits, result_mat)
def add(lhs: QOperator, rhs: QOperator):
q1 = lhs.qubits
q2 = rhs.qubits
qmax = np.max(np.concatenate((q1, q2))) + 1
lhs_big = lhs.broadcast(qmax)
rhs_big = rhs.broadcast(qmax)
result_mat = lhs_big + rhs_big
result_qubits = [i for i in range(qmax)]
return QOperator(result_qubits, result_mat)
# the identity operator I
def identity(qubits=[0]):
return QOperator(qubits, np.identity(2 ** len(qubits)))
def pauli_0(q=0):
return QOperator([q], np.identity(2))
# pauli's σ operators
def pauli_x(q = 0):
return QOperator([q], np.array([[0, 1], [1, 0]]))
def pauli_y(q = 0):
return QOperator([q], np.array([[0, 0 - 1j], [0 + 1j, 0]]))
def pauli_z(q = 0):
return QOperator([q], np.array([[1, 0], [0, -1]]))
def hadamard(q = 0):
return QOperator([q], 2**(-0.5) * np.array([[1, 1], [1, -1]], dtype=np.complex64))
σ_x = pauli_x
σ_y = pauli_y
σ_z = pauli_z
# which pauli operator?
def pauli(i, q = 0):
if i == 0:
return pauli_0(q)
elif i == 1:
return pauli_x(q)
elif i == 2:
return pauli_y(q)
elif i == 3:
return pauli_z(q)
else:
# should not happen
return pauli_0(q)
σ = pauli
# CNOT gate, 0 as control, 1 as target
def cnot(qubits = [0, 1]):
return QOperator(qubits, np.array([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0],
]))
# C-Z gate
def cphase(qubits = [0, 1]):
return QOperator(qubits, np.array([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, -1],
]))
def cpauli(pauli, qubits = [0, 1]):
if pauli == 'x':
return cnot(qubits)
elif pauli == 'z':
return cphase(qubits)
else:
# WARN this should not happen
return cnot(qubits)
import numpy as np
from functools import reduce
from utils import *
from QOperator import *
from QChannel import *
from DensityOperator import *
'''
q1, q2, q3, q4 are the 4 ancilla qubits
q2 ---------∎--⊤--M(x)
q1 ---∎--⊤--|--Z--M(x)
data1 ---|--σ--|---------
| |
q4 ---|-----∎--⊤--M(x)
q3 ---∎--⊤-----Z--M(x)
data2 ------σ------------
'''
def nickerson_1(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1, q2, q3 ,q4, pauli):
ρ1 = bell_pair(err_model.p_n, [q1, q3])
ρ2 = bell_pair(err_model.p_n, [q2, q4])
ρ.merge(ρ1)
ρ.merge(ρ2)
# 2 control-pauli gate
c1 = cpauli(pauli, [q1, data1])
c2 = cpauli(pauli, [q2, data2])
c3 = cphase([q2, q1])
c4 = cphase([q4, q3])
ρ.evolution(c1, err_model.p_g)
ρ.evolution(c2, err_model.p_g)
ρ.evolution(c3, err_model.p_g)
ρ.evolution(c4, err_model.p_g)
# measurement
ρ.bell_meausure(q1, q3, 'x', err_model.p_m)
ρ.bell_meausure(q1, q4, 'x', err_model.p_m)
return ρ
'''
q1, q2, q3, q4 are the 4 ancilla qubits
`stringent_plus` indicates whether the blank-separated region is performed
q2 ------∎--⊤--M(x)--∎--⊤--M(x)---- --∎--⊤--M(x)-- ------
q1 ---∎--|--X--------|--Z-------⊤-- --|--Z-------- --M(x)
data1 ---|--|-----------|----------σ-- --|----------- --------
| | | |
q4 ---|--∎--⊤--M(x)--∎--⊤--M(x)---- --∎--⊤--M(x)-- ------
q3 ---∎-----X-----------Z-------⊤-- -----Z-------- --M(x)
data2 -----------------------------σ-- -------------- --------
'''
def nickerson_2(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1, q2, q3 ,q4, pauli, stringent_plus = True):
ρ1 = bell_pair(err_model.p_n, [q1, q3])
ρ2 = bell_pair(err_model.p_n, [q2, q4])
ρ.merge(ρ1)
ρ.merge(ρ2)
c1 = cnot([q2, q1])
c2 = cnot([q4, q3])
ρ.evolution(c1, err_model.p_g)
ρ.evolution(c2, err_model.p_g)
ρ.bell_meausure(q2, q4, 'x', err_model.p_m)
ρ3 = bell_pair(err_model.p_n, [q2, q4])
ρ.merge(ρ2)
c3 = cphase([q2, q1])
c4 = cphase([q4, q3])
ρ.evolution(c3, err_model.p_g)
ρ.evolution(c4, err_model.p_g)
ρ.bell_meausure(q2, q4, 'x', err_model.p_m)
c5 = cpauli(pauli, [q1, data1])
c6 = cpauli(pauli, [q3, data2])
ρ.evolution(c5, err_model.p_g)
ρ.evolution(c6, err_model.p_g)
if stringent_plus:
ρ4 = bell_pair(err_model.p_n, [q2, q4])
ρ.merge(ρ4)
c7 = cphase([q2, q1])
c8 = cphase([q4, q3])
ρ.evolution(c7, err_model.p_g)
ρ.evolution(c8, err_model.p_g)
ρ.bell_meausure(q2, q4, 'x', err_model.p_m)
pass
ρ.bell_meausure(q1, q3, 'x', err_model.p_m)
'''
2 5
1 4
0 3
'''
def make_bell(err_model: ErrorModel, stringent = True, stringent_plus = True):
ρ = bell_pair(err_model.p_n, [0, 3])
nickerson_1(ρ, err_model, 0, 3, 1, 2, 4, 5, 'z')
nickerson_1(ρ, err_model, 0, 3, 1, 2, 4, 5, 'x')
if stringent:
nickerson_2(ρ, err_model, 0, 3, 1, 2, 4, 5, 'z', stringent_plus)
nickerson_2(ρ, err_model, 0, 3, 1, 2, 4, 5, 'x', stringent_plus)
return ρ
def make_GHZ(ρ1: DensityOperator, ρ2: DensityOperator, err_model: ErrorModel, stringent = True, stringent_plus = True):
pass
def stablizer_measurement(err_model: ErrorModel, stringent = True, stringent_plus = True):
ρ1 = make_bell(err_model, stringent, stringent_plus)
ρ2 = make_bell(err_model, stringent, stringent_plus)
ρ2.alter_qubits([6, 9])
import calc_channel as cc
import numpy as np
ρ = cc.DensityOperator(1, np.array([[1, 0], [0, 0]]))
channel = cc.depolarizing_channel(0.01)
ρ.channel(channel)
cnot = cc.cnot()
x = cc.σ_x()
z = cc.σ_z()
r = cc.multiply(z, cnot)
r = cc.multiply(r, z)
# print(ρ.operator)
print(r.operator)
'''
A: 0, 1, 2
B: 3, 4, 5
C: 6, 7, 8
D: 9, 10, 11
'''
import numpy as np
# get_bin_digits(6) -> array([0, 1, 1])
# get_bin_digits(6, 5) -> array([0, 1, 1, 0, 0])
def get_bin_digits(num, min_digits = None):
digits = np.array([], dtype=int)
while num >= 1:
digits = np.append(digits, num % 2)
num = num // 2
# to satisfy minimum digits
if min_digits is not None and len(digits) < min_digits:
digits = np.concatenate((digits, np.zeros(min_digits - len(digits), dtype=int)))
return digits
def get_n_digits(num, n, min_digits = None):
digits = np.array([], dtype=int)
while num >= 1:
digits = np.append(digits, num % n)
num = num // n
# to satisfy minimum digits
if min_digits is not None and len(digits) < min_digits:
digits = np.concatenate((digits, np.zeros(min_digits - len(digits), dtype=int)))
return digits
# [0, 1, 1] -> 6
def from_bin_digits(digits):
result = 0
for (i, d) in enumerate(digits):
result = result + d * 2 ** i
return result
def from_n_digits(digits, n):
result = 0
for (i, d) in enumerate(digits):
result = result + d * n ** i
return result
def contains(l, target):
for item in l:
if target == item:
return True
return False
class ErrorModel(object):
def __init__(self, p_n = 0, p_g = 0, p_m = 0):
self.p_n = p_n
self.p_g = p_g
self.p_m = p_m
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment