surface17_solution.py

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt

# We need to define for logical operators and for stabilizers the set qubits on
# which they act.
# This can be done using binary representation of an integer, in which the last
# bit indicates the physical qubit D1, and so on. For example, a logical
# X operator that acts on qubits D1, D4, and D7 is defined as:
xL = 0b001001001

# X stabilizer generators in symplectic notation
# qubit:987654321
s1 =  0b000000110
s2 =  0b000011011
s3 =  0b110110000
s4 =  0b011000000
stabilizerGenerators = [s1, s2, s3, s4]

def getCommutator(a, b):
	"""
	Check whether two operators made of Pauli operators commute or anticommute.
	a consists of Pauli X operators on a subset of the data qubits, b consists
	of Pauli Z operators on a different subset of data qubits. The data qubits
	are defined by the bits in the integers a and b.
	E.g. a=0b0101 has Pauli X operators on data qubits 1 and 3.
	Return 0 if a and b commute and 1 otherwise.
	"""
	# Because single qubit Paulis anticommute we just sum up (mod 2) the mutual
	# ones in the bitstings that represent the multiqubit operators in
	# symplectic notation

	# In python 3.10 you can uncomment the following line to make the code run faster
	#return (a & b).bit_count() % 2
	return bin(a & b)[2:].count('1')%2

def getSyndrome(configuration):
	"""
	Given a configuration of phase flip errors on the data qubits (defined by
	the binary representation of the integer "configuration"), compute the
	syndrome.
	Commutation/anticommutation with stabilizer n is marked as 0/1 in bit n
	of the syndrome. E.g. if an error anticommutes with stabilizers 1 and 3,
	the syndrome will be 0b0101.
	"""
	syndrome = 0
	for i, s in enumerate(stabilizerGenerators):
		syndrome += getCommutator(s, configuration)*2**i
	return syndrome

def getCorrection(syndrome):
	"""
	Given a syndrome, return the correction. The nth bit of the correction
	is 1 if data qubit n should obtain a phase flip in the correction, and 0 if
	data qubit n should not be altered.
	"""
	# Look Up Table: apply minimal weight correction based on syndrome information
	correction = []
	bin_correction = 0
	if syndrome == 0: return 0
	# TASK: add the correct syndromes and corrections
	if syndrome == 0b0001: correction = [3]
	if syndrome == 0b0010: correction = [1]
	if syndrome == 0b0011: correction = [2]
	if syndrome == 0b0100: correction = [9]
	if syndrome == 0b0101: correction = [3,9]
	if syndrome == 0b0110: correction = [5]
	if syndrome == 0b0111: correction = [2,9]
	if syndrome == 0b1000: correction = [7]
	if syndrome == 0b1001: correction = [3,7]
	if syndrome == 0b1010: correction = [1,7]
	if syndrome == 0b1011: correction = [2,7]
	if syndrome == 0b1100: correction = [8]
	if syndrome == 0b1101: correction = [3,8]
	if syndrome == 0b1110: correction = [1,8]
	if syndrome == 0b1111: correction = [2,8]
	for i in correction:
		bin_correction += 2**(i-1)
	return bin_correction

# set range of physical error rates to montecarlo sample
pRange = np.logspace(-3, 0, 7)

fail_log = []
logical_failure_rate = []
for p in pRange:
	fail_log.append([])
	# sample a reasonable amount of times given the probability of errors
	for _ in range(int(1000/p)):

		# start out with clean state and randomly place Z-errors on data qubits
		configuration = 0
		for i in range(9):
			if np.random.random() < p:
				configuration += 2**i

		# TASK: add syndrome readout and apply correction operators according to look up table
		syndrome = getSyndrome(configuration)
		correction = getCorrection(syndrome)
		after_correction = configuration ^ correction

		# TASK: check if the resulting configuration has caused a logical Z error
		if getCommutator(after_correction, xL) == 0:
			fail_log[-1].append(0)
		else:
			fail_log[-1].append(1)

		#print('configuration', format(configuration,'09b'))
		#print('syndrome', format(syndrome,'04b'))
		#print('correction', format(correction,'09b'))
		#print('result', format(after_correction,'09b'))
		#print('logical X', getCommutator(after_correction, xL))

	# calculate logical failure rate from samples
	logical_failure_rate.append(np.sum(fail_log[-1])/len(fail_log[-1]))

# Plot the results
plt.grid()
plt.xlabel('$p$')
plt.ylabel('$p_L$')
plt.title('Pseudothreshold of Surface-17')
plt.loglog(pRange, logical_failure_rate, '-x')
plt.loglog(pRange, pRange, 'k:', alpha=0.5)
plt.show()