From 69579745cae2eee38e03e66e01c3d8c9f9a44902 Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Fri, 24 Jun 2022 11:23:33 +0200
Subject: [PATCH] surface 17 depolarizing: adapted for SS22

---
 exercises/surface17_depolarizing.py           | 123 ++++++------
 exercises/surface17_depolarizing_optimized.py | 177 +++++++++---------
 exercises/surface17_optimized_exercise.py     |   8 +-
 exercises/surface17_optimized_solution.py     |   8 +-
 4 files changed, 156 insertions(+), 160 deletions(-)

diff --git a/exercises/surface17_depolarizing.py b/exercises/surface17_depolarizing.py
index 1ce7366..2c874bd 100644
--- a/exercises/surface17_depolarizing.py
+++ b/exercises/surface17_depolarizing.py
@@ -1,30 +1,37 @@
 #!/usr/bin/env python3
-'''
-Solution for exercise 13.1f, Monte Carlo simulation of surface-17
+"""
+Solution for exercise 12.3a, Monte Carlo simulation of surface-17
 quantum error correction for depolarizing noise.
 This script is much more efficient if the numba package is available.
-'''
+"""
 import numpy as np
+from random import random
 import matplotlib.pyplot as plt
 
+try:
+    from numba import jit, prange
+except ImportError:
+    prange = range
+    jit = lambda *args, **kwargs: (lambda x: x)
+
 # logical X and Z in symplectic notation
-xL = 0b001010100
-zL = 0b100010001
+xL = 0b001001001
+zL = 0b000000111
 # X stabilizer generators in symplectic notation
-sx1 = 0b000011011
-sx2 = 0b000100100
-sx3 = 0b001001000
-sx4 = 0b110110000
+sx1 = 0b000000110
+sx2 = 0b000011011
+sx3 = 0b110110000
+sx4 = 0b011000000
 # Z stabilizer generators in symplectic notation
-sz1 = 0b000000011
-sz2 = 0b000110110
-sz3 = 0b011011000
-sz4 = 0b110000000
-xstabilizerGenerators = (sx1,sx2,sx3,sx4)
-zstabilizerGenerators = (sz1,sz2,sz3,sz4)
+sz1 = 0b000001001
+sz2 = 0b011011000
+sz3 = 0b000110110
+sz4 = 0b100100000
+xstabilizerGenerators = (sx1, sx2, sx3, sx4)
+zstabilizerGenerators = (sz1, sz2, sz3, sz4)
 
 # Look-up table mapping 9 bit integers to the parity of the number of bits:
-# Definition: parity[x] = 1 if bitstring(x) has odd number of '1's, otherwise 0
+# Definition: parity[x] = 1 if bitstring(x) has odd number of '1's, otherwise 0.
 parity = bytes(bin(x).count('1') % 2 for x in range(2**9))
 
 # Z corrections detected by X stabilizers
@@ -32,21 +39,21 @@ zcorrections = (
         # binary data show which data qubit should be corrected:
         # 987654321 ← data qubit number
         0b000000000, # syndrom: -,       correction: -
-        0b000000001, # syndrom: 1,       correction: 1 (could also be 2)
-        0b000000100, # syndrom: 2,       correction: 3
-        0b000000101, # syndrom: 1,2,     correction: 1,3 (could also be 2,3 or 5,6)
-        0b001000000, # syndrom: 3,       correction: 7
-        0b000001000, # syndrom: 1,3,     correction: 4
-        0b001000100, # syndrom: 2,3,     correction: 3,7
-        0b000001100, # syndrom: 1,2,3,   correction: 3,4
-        0b100000000, # syndrom: 4,       correction: 9 (could also be 8)
-        0b000010000, # syndrom: 1,4,     correction: 5
-        0b000100000, # syndrom: 2,4,     correction: 6
-        0b000010100, # syndrom: 1,2,4,   correction: 3,5 (could also be 1,6)
-        0b101000000, # syndrom: 3,4,     correction: 7,9 (could also be 7,8 or 4,5)
-        0b001010000, # syndrom: 1,3,4,   correction: 5,7 (could also be 4,9 or 4,8)
-        0b001100000, # syndrom: 2,3,4,   correction: 6,7
-        0b000101000, # syndrom: 1,2,3,4, correction: 4,6
+        0b000000100, # syndrom: 1,       correction: 3
+        0b000000001, # syndrom: 2,       correction: 1 (could also be 4)
+        0b000000010, # syndrom: 1,2,     correction: 2
+        0b100000000, # syndrom: 3,       correction: 9 (could also be 6)
+        0b100000100, # syndrom: 1,3,     correction: 3,9 (could also be 3,6)
+        0b000010000, # syndrom: 2,3,     correction: 5
+        0b100000010, # syndrom: 1,2,3,   correction: 2,9 (could also be 2,6 or 3,5)
+        0b001000000, # syndrom: 4,       correction: 7
+        0b001000100, # syndrom: 1,4,     correction: 3,7
+        0b001000001, # syndrom: 2,4,     correction: 1,7 (could also be 4,7)
+        0b001000010, # syndrom: 1,2,4,   correction: 2,7
+        0b010000000, # syndrom: 3,4,     correction: 8
+        0b010000100, # syndrom: 1,3,4,   correction: 3,8
+        0b010000001, # syndrom: 2,3,4,   correction: 1,8 (could also be 4,8 or 5,7)
+        0b010000010, # syndrom: 1,2,3,4, correction: 2,8
         )
 
 # X corrections detected by Z stabilizers
@@ -55,60 +62,40 @@ xcorrections = (
         # 987654321 ← data qubit number
         0b000000000, # syndrom: -,       correction: -
         0b000000001, # syndrom: 1,       correction: 1
-        0b000000100, # syndrom: 2,       correction: 3 (could also be 6)
-        0b000000010, # syndrom: 1,2,     correction: 2
-        0b001000000, # syndrom: 3,       correction: 7 (could also be 4)
-        0b001000001, # syndrom: 1,3,     correction: 1,7 (could also be 1,4 or 2,5)
+        0b001000000, # syndrom: 2,       correction: 7 (could also be 8)
+        0b000001000, # syndrom: 1,2,     correction: 4
+        0b000000100, # syndrom: 3,       correction: 3 (could also be 2)
+        0b000000101, # syndrom: 1,3,     correction: 1,3 (could also be 1,2)
         0b000010000, # syndrom: 2,3,     correction: 5
-        0b000010001, # syndrom: 1,2,3,   correction: 1,5 (could also be 2,4 or 2,7)
+        0b000001100, # syndrom: 1,2,3,   correction: 3,4 (could also be 2,4 or 1,5)
         0b100000000, # syndrom: 4,       correction: 9
         0b100000001, # syndrom: 1,4,     correction: 1,9
-        0b100000100, # syndrom: 2,4,     correction: 3,9 (could also be 6,9 or 5,8)
-        0b100000010, # syndrom: 1,2,4,   correction: 2,9
-        0b010000000, # syndrom: 3,4,     correction: 8
-        0b010000001, # syndrom: 1,3,4,   correction: 1,8
-        0b100010000, # syndrom: 2,3,4,   correction: 5,9 (could also be 3,8 or 6,8)
-        0b010000010, # syndrom: 1,2,3,4, correction: 2,8
+        0b101000000, # syndrom: 2,4,     correction: 7,9 (could also be 8,9)
+        0b100001000, # syndrom: 1,2,4,   correction: 4,9
+        0b000100000, # syndrom: 3,4,     correction: 6
+        0b000100001, # syndrom: 1,3,4,   correction: 1,6
+        0b001100000, # syndrom: 2,3,4,   correction: 6,7 (could also be 6,8 or 5,9)
+        0b000101000, # syndrom: 1,2,3,4, correction: 4,6
         )
 
-
-################################################################
-# THIS BLOCK COULD BE OMITTED, but it can significantly improve the performance
-# if the numba package is available. Numba provides just-in-time compilation
-# and parallelization that speed up the for loop in montecarlo().
-try:
-    # Try to import just-in-time compilation function decorator from numba.
-    from numba import jit, prange
-    NUMBA = True
-except ImportError:
-    # prange behaves like range, but allows for parallelization. Replace it
-    # with range if numba is not available.
-    prange = range
-    # If numba is not available, jit should take a dummy value.
-    jit = lambda *args, **kwargs: (lambda x: x)
-    NUMBA = False
-# The decorator @jit causes this function to be compiled and not just
-# parsed by the python interpreter:
 @jit(nopython=True, parallel=True, cache=False)
-################################################################
 def montecarlo(p, samples):
-    '''
+    """
     Simulate surface 17 code with error rate p/3 for X, Y and Z errors on data
     qubits for given number of samples. Return number of logical X, Y, and Z
     errors.
-    '''
+    """
     xfailures = 0
     yfailures = 0
     zfailures = 0
     # sample a reasonable amount of times given the probability of errors
     for _ in prange(samples):
-
         # start out with clean state and randomly place X, Z and XZ errors on data qubits
         xerrors = 0
         zerrors = 0
         # Iterate over the 9 qubits
         for i in range(9):
-            r = np.random.random()
+            r = random()
             if r < p/3:
                 xerrors |= 2**i
             elif r < 2*p/3:
@@ -144,10 +131,10 @@ def montecarlo(p, samples):
 
 
 def main(relative_samples=1000, resolution=16):
-    '''
+    """
     Simulate and plot logical error rates of surface 17 code with depolarizing
     noise.
-    '''
+    """
     # set range of physical error rates to MC sample
     pRange = np.logspace(-3, 0, resolution)
     # define number of samples for each physical error rate
@@ -184,4 +171,4 @@ def main(relative_samples=1000, resolution=16):
 
 
 if __name__ == '__main__':
-    main(10000 if NUMBA else 1000)
+    main(1000)
diff --git a/exercises/surface17_depolarizing_optimized.py b/exercises/surface17_depolarizing_optimized.py
index ebc0b7b..b713142 100644
--- a/exercises/surface17_depolarizing_optimized.py
+++ b/exercises/surface17_depolarizing_optimized.py
@@ -1,120 +1,121 @@
 #!/usr/bin/env python3
-'''
-Solution for exercise 13.1f, Monte Carlo simulation of surface-17
-quantum error correction for depolarizing noise with some optimizations.
-This script is much more efficient if the numba package is available.
-'''
+"""
+Solution for exercise 12.3a, Monte Carlo simulation of surface-17
+quantum error correction for depolarizing noise.
+"""
 import numpy as np
-from random import random
 import matplotlib.pyplot as plt
+from random import random
+from numba import jit, prange
 
 # logical X and Z in symplectic notation
-xL = 0b001010100
-zL = 0b100010001
+xL = 0b001001001
+zL = 0b000000111
 # X stabilizer generators in symplectic notation
-sx1 = 0b000011011
-sx2 = 0b000100100
-sx3 = 0b001001000
-sx4 = 0b110110000
+sx1 = 0b000000110
+sx2 = 0b000011011
+sx3 = 0b110110000
+sx4 = 0b011000000
 # Z stabilizer generators in symplectic notation
-sz1 = 0b000000011
-sz2 = 0b000110110
-sz3 = 0b011011000
-sz4 = 0b110000000
+sz1 = 0b000001001
+sz2 = 0b011011000
+sz3 = 0b000110110
+sz4 = 0b100100000
 xstabilizerGenerators = (sx1,sx2,sx3,sx4)
 zstabilizerGenerators = (sz1,sz2,sz3,sz4)
 
-parity = bytes(bin(x).count('1') % 2 for x in range(1<<9))
+# Look-up table mapping 9 bit integers to the parity of the number of bits:
+# Definition: parity[x] = 1 if bitstring(x) has odd number of '1's, otherwise 0
+parity = bytes(x.bit_count() % 2 for x in range(1<<9))
 
 # Z corrections detected by X stabilizers
 zcorrections = (
+        # binary data show which data qubit should be corrected:
+        # 987654321 ← data qubit number
         0b000000000, # syndrom: -,       correction: -
-        0b000000001, # syndrom: 1,       correction: 1 (could also be 2)
-        0b000000100, # syndrom: 2,       correction: 3
-        0b000000101, # syndrom: 1,2,     correction: 1,3 (could also be 2,3 or 5,6)
-        0b001000000, # syndrom: 3,       correction: 7
-        0b000001000, # syndrom: 1,3,     correction: 4
-        0b001000100, # syndrom: 2,3,     correction: 3,7
-        0b000001100, # syndrom: 1,2,3,   correction: 3,4
-        0b100000000, # syndrom: 4,       correction: 9 (could also be 8)
-        0b000010000, # syndrom: 1,4,     correction: 5
-        0b000100000, # syndrom: 2,4,     correction: 6
-        0b000010100, # syndrom: 1,2,4,   correction: 3,5 (could also be 1,6)
-        0b101000000, # syndrom: 3,4,     correction: 7,9 (could also be 7,8 or 4,5)
-        0b001010000, # syndrom: 1,3,4,   correction: 5,7 (could also be 4,9 or 4,8)
-        0b001100000, # syndrom: 2,3,4,   correction: 6,7
-        0b000101000, # syndrom: 1,2,3,4, correction: 4,6
+        0b000000100, # syndrom: 1,       correction: 3
+        0b000000001, # syndrom: 2,       correction: 1 (could also be 4)
+        0b000000010, # syndrom: 1,2,     correction: 2
+        0b100000000, # syndrom: 3,       correction: 9 (could also be 6)
+        0b100000100, # syndrom: 1,3,     correction: 3,9 (could also be 3,6)
+        0b000010000, # syndrom: 2,3,     correction: 5
+        0b100000010, # syndrom: 1,2,3,   correction: 2,9 (could also be 2,6 or 3,5)
+        0b001000000, # syndrom: 4,       correction: 7
+        0b001000100, # syndrom: 1,4,     correction: 3,7
+        0b001000001, # syndrom: 2,4,     correction: 1,7 (could also be 4,7)
+        0b001000010, # syndrom: 1,2,4,   correction: 2,7
+        0b010000000, # syndrom: 3,4,     correction: 8
+        0b010000100, # syndrom: 1,3,4,   correction: 3,8
+        0b010000001, # syndrom: 2,3,4,   correction: 1,8 (could also be 4,8 or 5,7)
+        0b010000010, # syndrom: 1,2,3,4, correction: 2,8
         )
 
 # X corrections detected by Z stabilizers
 xcorrections = (
+        # binary data show which data qubit should be corrected:
+        # 987654321 ← data qubit number
         0b000000000, # syndrom: -,       correction: -
         0b000000001, # syndrom: 1,       correction: 1
-        0b000000100, # syndrom: 2,       correction: 3 (could also be 6)
-        0b000000010, # syndrom: 1,2,     correction: 2
-        0b001000000, # syndrom: 3,       correction: 7 (could also be 4)
-        0b001000001, # syndrom: 1,3,     correction: 1,7 (could also be 1,4 or 2,5)
+        0b001000000, # syndrom: 2,       correction: 7 (could also be 8)
+        0b000001000, # syndrom: 1,2,     correction: 4
+        0b000000100, # syndrom: 3,       correction: 3 (could also be 2)
+        0b000000101, # syndrom: 1,3,     correction: 1,3 (could also be 1,2)
         0b000010000, # syndrom: 2,3,     correction: 5
-        0b000010001, # syndrom: 1,2,3,   correction: 1,5 (could also be 2,4 or 2,7)
+        0b000001100, # syndrom: 1,2,3,   correction: 3,4 (could also be 2,4 or 1,5)
         0b100000000, # syndrom: 4,       correction: 9
         0b100000001, # syndrom: 1,4,     correction: 1,9
-        0b100000100, # syndrom: 2,4,     correction: 3,9 (could also be 6,9 or 5,8)
-        0b100000010, # syndrom: 1,2,4,   correction: 2,9
-        0b010000000, # syndrom: 3,4,     correction: 8
-        0b010000001, # syndrom: 1,3,4,   correction: 1,8
-        0b100010000, # syndrom: 2,3,4,   correction: 5,9 (could also be 3,8 or 6,8)
-        0b010000010, # syndrom: 1,2,3,4, correction: 2,8
+        0b101000000, # syndrom: 2,4,     correction: 7,9 (could also be 8,9)
+        0b100001000, # syndrom: 1,2,4,   correction: 4,9
+        0b000100000, # syndrom: 3,4,     correction: 6
+        0b000100001, # syndrom: 1,3,4,   correction: 1,6
+        0b001100000, # syndrom: 2,3,4,   correction: 6,7 (could also be 6,8 or 5,9)
+        0b000101000, # syndrom: 1,2,3,4, correction: 4,6
         )
 
-direct_xcorrections = bytes(parity[(c ^ xcorrections[sum(parity[s&c]<<i for i,s in enumerate(zstabilizerGenerators))]) & zL] for c in range(1<<9))
-direct_zcorrections = bytes(parity[(c ^ zcorrections[sum(parity[s&c]<<i for i,s in enumerate(xstabilizerGenerators))]) & xL] for c in range(1<<9))
+direct_xcorrections = bytearray(1<<9)
+for configuration in range(1<<9):
+    syndrome = 0
+    for i,s in enumerate(zstabilizerGenerators):
+        syndrome |= parity[s & configuration] << i
+    corrected = configuration ^ xcorrections[syndrome]
+    direct_xcorrections[configuration] = parity[corrected & zL]
+direct_xcorrections = bytes(direct_xcorrections)
+
+direct_zcorrections = bytearray(1<<9)
+for configuration in range(1<<9):
+    syndrome = 0
+    for i,s in enumerate(xstabilizerGenerators):
+        syndrome |= parity[s & configuration] << i
+    corrected = configuration ^ zcorrections[syndrome]
+    direct_zcorrections[configuration] = parity[corrected & xL]
+direct_zcorrections = bytes(direct_zcorrections)
 
-################################################################
-# THIS BLOCK COULD BE OMITTED, but it can significantly improve the performance
-# if the numba package is available. Numba provides just-in-time compilation
-# and parallelization that speed up the for loop in montecarlo().
-try:
-    # Try to import just-in-time compilation function decorator from numba.
-    from numba import jit, prange
-except ImportError:
-    # prange behaves like range, but allows for parallelization. Replace it
-    # with range if numba is not available.
-    prange = range
-    # If numba is not available, jit should take a dummy value.
-    jit = lambda *args, **kwargs: (lambda x: x)
-# The decorator @jit causes this function to be compiled and not just
-# parsed by the python interpreter:
 @jit(nopython=True, parallel=True, cache=False)
-################################################################
 def montecarlo(p, samples):
-    '''
+    """
     Simulate surface 17 code with error rate p/3 for X, Y and Z errors on data
     qubits for given number of samples. Return number of logical X, Y, and Z
     errors.
-    '''
+    """
     xfailures = 0
     yfailures = 0
     zfailures = 0
     # sample a reasonable amount of times given the probability of errors
     for _ in prange(samples):
-
-        # start out with clean state and randomly place Z-errors on data qubits
+        # start out with clean state and randomly place X, Z and XZ errors on data qubits
         xerrors = 0
         zerrors = 0
+        # Iterate over the 9 qubits
         for i in range(9):
             r = random()
-            #if r < p:
-            #    zerrors |= 1 << i
             if r < p/3:
-                xerrors |= 1 << i
+                xerrors |= 1<<i
             elif r < 2*p/3:
-                zerrors |= 1 << i
+                zerrors |= 1<<i
             elif r < p:
-                xerrors |= 1 << i
-                zerrors |= 1 << i
+                xerrors |= 1<<i
+                zerrors |= 1<<i
 
-        # check if the resulting configuration has caused a logical Z error
-        # logical Z should always anticommute with logical X
         x_logical_error = direct_xcorrections[xerrors]
         z_logical_error = direct_zcorrections[zerrors]
         if x_logical_error and z_logical_error:
@@ -126,21 +127,29 @@ def montecarlo(p, samples):
     return xfailures, yfailures, zfailures
 
 
-def main(relative_samples=10000, resolution=16):
-    '''
+def main(relative_samples=1000, resolution=16):
+    """
     Simulate and plot logical error rates of surface 17 code with depolarizing
     noise.
-    '''
+    """
     # set range of physical error rates to MC sample
     pRange = np.logspace(-3, 0, resolution)
     # define number of samples for each physical error rate
     sampleRange = np.int_(relative_samples/pRange)
 
+    x_failure_rates = np.empty_like(pRange)
+    y_failure_rates = np.empty_like(pRange)
+    z_failure_rates = np.empty_like(pRange)
     # run the MC simulation
-    logical_error_counts = np.array([montecarlo(*ps) for ps in zip(pRange, sampleRange)])
-    # calculate logical failure rate from samples
-    logical_error_rates = logical_error_counts / sampleRange.reshape((resolution,1))
+    for i in range(resolution):
+        xfailures, yfailures, zfailures = montecarlo(pRange[i], sampleRange[i])
+        x_failure_rates[i] = xfailures / sampleRange[i]
+        y_failure_rates[i] = yfailures / sampleRange[i]
+        z_failure_rates[i] = zfailures / sampleRange[i]
+
+    total_failure_rates = x_failure_rates + y_failure_rates + z_failure_rates
 
+    # Plot the results
     fig, ax = plt.subplots()
     ax.grid()
     ax.set_xlabel('$p$')
@@ -148,10 +157,10 @@ def main(relative_samples=10000, resolution=16):
     ax.set_title('Pseudothreshold of Surface-17, depolarizing noise')
     ax.set_xscale('log')
     ax.set_yscale('log')
-    ax.plot(pRange, logical_error_rates[:,0], '-x', label='logical X error')
-    ax.plot(pRange, logical_error_rates[:,1], '-x', label='logical Y error')
-    ax.plot(pRange, logical_error_rates[:,2], '-x', label='logical Z error')
-    ax.plot(pRange, logical_error_rates.sum(axis=1), '-x', label='any logical error')
+    ax.plot(pRange, x_failure_rates, '-x', label='logical X error')
+    ax.plot(pRange, y_failure_rates, '-x', label='logical Y error')
+    ax.plot(pRange, z_failure_rates, '-x', label='logical Z error')
+    ax.plot(pRange, total_failure_rates, '-x', label='any logical error')
     ax.plot(pRange, pRange, 'k:', alpha=0.5)
     ax.legend()
     fig.tight_layout()
@@ -159,4 +168,4 @@ def main(relative_samples=10000, resolution=16):
 
 
 if __name__ == '__main__':
-    main()
+    main(10000)
diff --git a/exercises/surface17_optimized_exercise.py b/exercises/surface17_optimized_exercise.py
index 12d92bf..5a94660 100644
--- a/exercises/surface17_optimized_exercise.py
+++ b/exercises/surface17_optimized_exercise.py
@@ -63,14 +63,14 @@ correction_table = (
 
 
 @jit(nopython=True, cache=False, parallel=True)
-def simulate(p, runs):
+def simulate(p, samples):
     """
-    Simulate n runs with error probability p.
+    Simulate n samples 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):
+    for _ in prange(samples):
 
         # start out with clean state and randomly place Z-errors on data qubits
         configuration = 0
@@ -88,7 +88,7 @@ def simulate(p, runs):
             failures += 1
 
     # calculate logical failure rate from samples
-    return failures / runs
+    return failures / samples
 
 def main():
     # set range of physical error rates to MC sample
diff --git a/exercises/surface17_optimized_solution.py b/exercises/surface17_optimized_solution.py
index 171b81f..f260daa 100644
--- a/exercises/surface17_optimized_solution.py
+++ b/exercises/surface17_optimized_solution.py
@@ -63,14 +63,14 @@ correction_table = (
 
 
 @jit(nopython=True, cache=False, parallel=True)
-def simulate(p, runs):
+def simulate(p, samples):
     """
-    Simulate n runs with error probability p.
+    Simulate n samples 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):
+    for _ in prange(samples):
 
         # start out with clean state and randomly place Z-errors on data qubits
         configuration = 0
@@ -88,7 +88,7 @@ def simulate(p, runs):
             failures += 1
 
     # calculate logical failure rate from samples
-    return failures / runs
+    return failures / samples
 
 def main():
     # set range of physical error rates to MC sample
-- 
GitLab