From 1d4ed796e0ac9d513e1ec0d201c4af8ea77cc81c Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Fri, 24 Jun 2022 11:10:12 +0200
Subject: [PATCH] surface 17 dephasing: adapted for SS22
---
exercises/surface17_exercise.py | 130 +++++++++++++---------
exercises/surface17_optimized_exercise.py | 111 ++++++++++++++++++
exercises/surface17_optimized_solution.py | 111 ++++++++++++++++++
exercises/surface17_solution.py | 130 ++++++++++++----------
4 files changed, 372 insertions(+), 110 deletions(-)
create mode 100644 exercises/surface17_optimized_exercise.py
create mode 100644 exercises/surface17_optimized_solution.py
diff --git a/exercises/surface17_exercise.py b/exercises/surface17_exercise.py
index bbbf8cc..723e7df 100644
--- a/exercises/surface17_exercise.py
+++ b/exercises/surface17_exercise.py
@@ -2,100 +2,124 @@
import numpy as np
import matplotlib.pyplot as plt
-# logical X in symplectic notation
-xL = 0b001010100
+# 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
-s1 = 0b000011011
-s2 = 0b000100100
-s3 = 0b001001000
-s4 = 0b110110000
-stabilizerGenerators = [s1,s2,s3,s4]
+# 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
-def getCommutator(a,b):
- # For two multiqubit Pauli X- or Z- operators we want to know
- # whether they commute (0) or anticommute (1). 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):
- # commutation/anticommutation with a stabilizer is marked as 0/1 in syndrome
+ """
+ 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 in range(len(stabilizerGenerators)):
- syndrome += getCommutator(stabilizerGenerators[i], configuration)*2**i
+ 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 == 0b0001: correction = [0]
- if syndrome == 0b1110: correction = [5,6]
- if syndrome == 0b0010: correction = [2]
- if syndrome == 0b1000: correction = [8]
- if syndrome == 0b0100: correction = [6]
+ if syndrome == 0: return 0
# TASK: add the correct syndromes and corrections
- if syndrome == ...: correction = [3]
- if syndrome == ...: correction = [2,4]
- if syndrome == 0b1101: correction = ...
+ 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 = ...
+ if syndrome == 0b0111: correction = [2,9]
+ if syndrome == 0b1000: correction = ...
if syndrome == 0b1001: correction = ...
- if syndrome == 0b1010: correction = [5]
- if syndrome == 0b1100: correction = [6,8]
- if syndrome == 0b0111: correction = [5,6]
- if syndrome == 0b0011: correction = [0,2]
- if syndrome == 0b0110: correction = [2,6]
- if syndrome == 0b1111: correction = [3,5]
+ if syndrome == 0b1010: correction = [1,7]
+ if syndrome == 0b1011: correction = [2,7]
+ if syndrome == 0b1100: correction = ...
+ 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
+ bin_correction += 2**(i-1)
return bin_correction
-# set range of physical error rates to MC sample
-pRange = np.logspace(-3,0,7)
+# set range of physical error rates to montecarlo sample
+pRange = np.logspace(-3, 0, 7)
fail_log = []
logical_failure_rate = []
-for p in range(len(pRange)):
-
+for p in pRange:
fail_log.append([])
# sample a reasonable amount of times given the probability of errors
- for _ in range(int(1000/pRange[p])):
-
+ 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() < pRange[p]: configuration += 2**i
-
+ if np.random.random() < p:
+ configuration += 2**i
+
# TASK: add syndrome readout and apply correction operators according to look up table
- syndrome = ...
- correction = ...
- after_correction = configuration^correction
+ syndrome = ...
+ correction = ...
+ after_correction = configuration ^ correction
# TASK: check if the resulting configuration has caused a logical Z error
if ... :
fail_log[-1].append(0)
else:
fail_log[-1].append(1)
-
+
+ # You can uncomment these lines to debug your code
#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]))
-#print('physical error rates', pRange)
-#print('failure', logical_failure_rate)
-
-#plt.rcParams['text.usetex'] = True
-#plt.rcParams.update({'font.size': 18})
+ # 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(np.logspace(-3,0,4),np.logspace(-3,0,4),'k:',alpha=0.5)
-plt.tight_layout()
+plt.loglog(pRange, logical_failure_rate, '-x')
+plt.loglog(pRange, pRange, 'k:', alpha=0.5)
plt.show()
diff --git a/exercises/surface17_optimized_exercise.py b/exercises/surface17_optimized_exercise.py
new file mode 100644
index 0000000..12d92bf
--- /dev/null
+++ b/exercises/surface17_optimized_exercise.py
@@ -0,0 +1,111 @@
+#!/usr/bin/env python3
+import numpy as np
+from random import random
+import matplotlib.pyplot as plt
+from numba import jit, prange
+
+# 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)
+
+# Lookup table of number of bits:
+# parity_table[x] = 1 if x has an odd number of ones, 0 otherwise
+# This implies getCommutator(a, b) == parity_table[a & b]
+parity_table = bytes(bin(x).count('1')%2 for x in range(1<<9))
+
+@jit(nopython=True)
+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 |= parity_table[s & configuration] << i
+ return syndrome
+
+# Lookup table of corrections:
+# correction_table[syndrome] == getCorrection(syndrome)
+correction_table = (
+ # TASK: add the correct syndromes and corrections
+ 0b000000000, # syndmome = 0, correction = []
+ 0b000000100, # syndrome = 0b0001, correction = [3]
+ 0b000000001, # syndrome = 0b0010, correction = [1]
+ 0b000000010, # syndrome = 0b0011, correction = [2]
+ 0b100000000, # syndrome = 0b0100, correction = [9]
+ 0b100000100, # syndrome = 0b0101, correction = [3,9]
+ ... , # syndrome = 0b0110, correction = ...
+ 0b100000010, # syndrome = 0b0111, correction = [2,9]
+ ... , # syndrome = 0b1000, correction = ...
+ ... , # syndrome = 0b1001, correction = ...
+ 0b001000001, # syndrome = 0b1010, correction = [1,7]
+ 0b001000010, # syndrome = 0b1011, correction = [2,7]
+ ... , # syndrome = 0b1100, correction = ...
+ 0b010000100, # syndrome = 0b1101, correction = [3,8]
+ 0b010000001, # syndrome = 0b1110, correction = [1,8]
+ 0b010000010, # syndrome = 0b1111, correction = [2,8]
+ )
+
+
+@jit(nopython=True, cache=False, parallel=True)
+def simulate(p, runs):
+ """
+ Simulate n runs with error probability p.
+ Return logical failure rate.
+ """
+ failures = 0
+ # sample a reasonable amount of times given the probability of errors
+ for _ in prange(runs):
+
+ # start out with clean state and randomly place Z-errors on data qubits
+ configuration = 0
+ for i in range(9):
+ if random() < p:
+ configuration |= 1<<i
+
+ # TASK: add syndrome readout and apply correction operators according to look up table
+ syndrome = ...
+ correction = ...
+ after_correction = configuration ^ correction
+
+ # TASK: check if the resulting configuration has caused a logical Z error
+ if ... :
+ failures += 1
+
+ # calculate logical failure rate from samples
+ return failures / runs
+
+def main():
+ # set range of physical error rates to MC sample
+ pRange = np.logspace(-3, 0, 21)
+
+ # Run the simulation
+ logical_failure_rate = [simulate(p, int(5000/p)) for p in pRange]
+
+ # 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()
+
+
+if __name__ == '__main__':
+ main()
diff --git a/exercises/surface17_optimized_solution.py b/exercises/surface17_optimized_solution.py
new file mode 100644
index 0000000..171b81f
--- /dev/null
+++ b/exercises/surface17_optimized_solution.py
@@ -0,0 +1,111 @@
+#!/usr/bin/env python3
+import numpy as np
+from random import random
+import matplotlib.pyplot as plt
+from numba import jit, prange
+
+# 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)
+
+# Lookup table of number of bits:
+# parity_table[x] = 1 if x has an odd number of ones, 0 otherwise
+# This implies getCommutator(a, b) == parity_table[a & b]
+parity_table = bytes(bin(x).count('1')%2 for x in range(1<<9))
+
+@jit(nopython=True)
+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 |= parity_table[s & configuration] << i
+ return syndrome
+
+# Lookup table of corrections:
+# correction_table[syndrome] == getCorrection(syndrome)
+correction_table = (
+ # TASK: add the correct syndromes and corrections
+ 0b000000000, # syndmome = 0, correction = []
+ 0b000000100, # syndrome = 0b0001, correction = [3]
+ 0b000000001, # syndrome = 0b0010, correction = [1]
+ 0b000000010, # syndrome = 0b0011, correction = [2]
+ 0b100000000, # syndrome = 0b0100, correction = [9]
+ 0b100000100, # syndrome = 0b0101, correction = [3,9]
+ 0b000010000, # syndrome = 0b0110, correction = [5]
+ 0b100000010, # syndrome = 0b0111, correction = [2,9]
+ 0b001000000, # syndrome = 0b1000, correction = [7]
+ 0b001000100, # syndrome = 0b1001, correction = [3,7]
+ 0b001000001, # syndrome = 0b1010, correction = [1,7]
+ 0b001000010, # syndrome = 0b1011, correction = [2,7]
+ 0b010000000, # syndrome = 0b1100, correction = [8]
+ 0b010000100, # syndrome = 0b1101, correction = [3,8]
+ 0b010000001, # syndrome = 0b1110, correction = [1,8]
+ 0b010000010, # syndrome = 0b1111, correction = [2,8]
+ )
+
+
+@jit(nopython=True, cache=False, parallel=True)
+def simulate(p, runs):
+ """
+ Simulate n runs with error probability p.
+ Return logical failure rate.
+ """
+ failures = 0
+ # sample a reasonable amount of times given the probability of errors
+ for _ in prange(runs):
+
+ # start out with clean state and randomly place Z-errors on data qubits
+ configuration = 0
+ for i in range(9):
+ if random() < p:
+ configuration |= 1<<i
+
+ # TASK: add syndrome readout and apply correction operators according to look up table
+ syndrome = getSyndrome(configuration)
+ correction = correction_table[syndrome]
+ after_correction = configuration ^ correction
+
+ # TASK: check if the resulting configuration has caused a logical Z error
+ if parity_table[after_correction & xL]:
+ failures += 1
+
+ # calculate logical failure rate from samples
+ return failures / runs
+
+def main():
+ # set range of physical error rates to MC sample
+ pRange = np.logspace(-3, 0, 21)
+
+ # Run the simulation
+ logical_failure_rate = [simulate(p, int(5000/p)) for p in pRange]
+
+ # 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()
+
+
+if __name__ == '__main__':
+ main()
diff --git a/exercises/surface17_solution.py b/exercises/surface17_solution.py
index fc3c5ef..a7f1527 100644
--- a/exercises/surface17_solution.py
+++ b/exercises/surface17_solution.py
@@ -1,83 +1,105 @@
#!/usr/bin/env python3
-'''
-Solution for exercise 13.1f, Monte Carlo simulation of surface-17 error
-quantum error correction for dephasing noise.
-A more efficient implementation of this can be found in the solution for
-depolarizing noise.
-'''
import numpy as np
import matplotlib.pyplot as plt
-# logical X in symplectic notation
-xL = 0b001010100
+# 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
-s1 = 0b000011011
-s2 = 0b000100100
-s3 = 0b001001000
-s4 = 0b110110000
-stabilizerGenerators = [s1,s2,s3,s4]
+# 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
-def getCommutator(a,b):
- # For two multiqubit Pauli X- or Z- operators we want to know
- # whether they commute (0) or anticommute (1). 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):
- # commutation/anticommutation with a stabilizer is marked as 0/1 in syndrome
+ """
+ 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 in range(len(stabilizerGenerators)):
- syndrome += getCommutator(stabilizerGenerators[i], configuration)*2**i
+ 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 == 0b0001: correction = [0] # could also be [1]
- if syndrome == 0b1110: correction = [5,6]
- if syndrome == 0b0010: correction = [2]
- if syndrome == 0b1000: correction = [8] # could also be [7]
- if syndrome == 0b0100: correction = [6]
- if syndrome == 0b0101: correction = [3]
- if syndrome == 0b1011: correction = [2,4]
- if syndrome == 0b1101: correction = [4,6] # could also be [3,8]
- if syndrome == 0b1001: correction = [4]
- if syndrome == 0b1010: correction = [5]
- if syndrome == 0b1100: correction = [6,8]
- if syndrome == 0b0111: correction = [5,6]
- if syndrome == 0b0011: correction = [0,2]
- if syndrome == 0b0110: correction = [2,6]
- if syndrome == 0b1111: correction = [3,5]
+ 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
+ bin_correction += 2**(i-1)
return bin_correction
-# set range of physical error rates to MC sample
-pRange = np.logspace(-3,0,7)
+# set range of physical error rates to montecarlo sample
+pRange = np.logspace(-3, 0, 7)
fail_log = []
logical_failure_rate = []
-for p in range(len(pRange)):
-
+for p in pRange:
fail_log.append([])
# sample a reasonable amount of times given the probability of errors
- for _ in range(int(1000/pRange[p])):
+ 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() < pRange[p]: configuration += 2**i
+ if np.random.random() < p:
+ configuration += 2**i
- # read syndrome and apply correction operators according to look up table
+ # TASK: add syndrome readout and apply correction operators according to look up table
syndrome = getSyndrome(configuration)
correction = getCorrection(syndrome)
- after_correction = configuration^correction
+ after_correction = configuration ^ correction
- # check if the resulting configuration has caused a logical Z error
- # we notice that we induced a logical Z by applying the correction
- # by checking if the resulting multiqubit Z-operator commutes with logical X
+ # TASK: check if the resulting configuration has caused a logical Z error
if getCommutator(after_correction, xL) == 0:
fail_log[-1].append(0)
else:
@@ -92,17 +114,11 @@ for p in range(len(pRange)):
# calculate logical failure rate from samples
logical_failure_rate.append(np.sum(fail_log[-1])/len(fail_log[-1]))
-#print('physical error rates', pRange)
-#print('failure', logical_failure_rate)
-
-#plt.rcParams['text.usetex'] = True
-#plt.rcParams.update({'font.size': 18})
-
+# 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(np.logspace(-3,0,4),np.logspace(-3,0,4),'k:',alpha=0.5)
-plt.tight_layout()
+plt.loglog(pRange, logical_failure_rate, '-x')
+plt.loglog(pRange, pRange, 'k:', alpha=0.5)
plt.show()
--
GitLab