LB: Added bell_plotting_helpers.py and s17_plotting_helpers.py, where we can...

LB: Added bell_plotting_helpers.py and s17_plotting_helpers.py, where we can share functionality for ploting. The notebooks are also pushed once, but bell_ploting.ipynb and s17_plotting.ipynb should be kept local (add to .gitignore)

plot/s17_plotting_helpers.py

+import stim
+import numpy as np
+import random
+import matplotlib.pyplot as plt
+import scipy
+
+
+
+# amplitude decay probability before measurement. 
+
+def simulate_circuit(circuit_in, 
+                     p_amp_ge = 0.01, 
+                     ancilla_fb = False, 
+                     after_clifford_depolarization_2 = 0.01): 
+    ### Operator string. 
+    # qb: 1, 3, 5
+    logical_qubits = [1,3,5]  # qubits of logical operator on S17
+    operator_string = ""
+    for i in range(circuit_in.num_qubits):
+        if i in logical_qubits: 
+            operator_string += "X"
+        else: 
+            operator_string += "I"
+            
+    # print(operator_string)
+
+    simulator = stim.TableauSimulator()
+
+    z_ancillas = [9,13,14,18]
+    x_ancillas = [2,16,11,25]
+    
+    z_ancillas_dict = {i: {'msmt_last': None } for i in z_ancillas}
+    x_ancillas_dict = {i: {'msmt_last': None } for i in x_ancillas}
+        
+    circuit = circuit_in.flattened() #circuit.flattened()  
+    for k in range(len(circuit)) :
+        instruction = circuit[k]
+        sub_circuit = circuit[k:k+1]
+        #print(instruction.name)
+        if  instruction.name == 'MR' :
+    
+            ts = instruction.targets_copy()
+            z_ancillas_list = []
+            x_ancillas_list = []
+            qubits = [q.value for q in ts]
+            
+            for qubit  in [q.value for q in  ts]:   
+                # apply flip: 
+                flip_applied = False
+                if ancilla_fb: 
+                    if qubit in z_ancillas:
+                        # if z_ancillas_dict[qubit]['msmt_last'] == [True for i in range(2)]: 
+                        # if z_ancillas_dict[qubit]['msmt_last'] == [True]: 
+                        if z_ancillas_dict[qubit]['msmt_last'] in [ [True, True, True], [False, True, True], [True, False, True]]: #
+                            simulator.x(qubit)
+                            flip_applied = True
+                    if qubit in x_ancillas:
+                        # if x_ancillas_dict[qubit]['msmt_last'] == [True for i in range(2)]: 
+                        if x_ancillas_dict[qubit]['msmt_last'] == [ [True, True, True], [False, True, True], [True, False, True]]:  #, [True, True, False]
+                        # if x_ancillas_dict[qubit]['msmt_last'] == [True]: 
+                            simulator.x(qubit) 
+                            # print('flip_applied',qubit, simulator.peek_bloch(qubit))
+                            flip_applied =True
+
+                ## apply noise    
+                    # check whether we are in |1>
+                if simulator.peek_bloch(qubit) == stim.PauliString("-Z"): 
+                    if random.random() < p_amp_ge: 
+                        simulator.x(qubit)
+                        
+                # we cannot change the measurement record... fake flip before msmt. 
+                # apply flip: 
+                if ancilla_fb and flip_applied: 
+                    simulator.x(qubit)
+                    flip_applied = False
+
+
+            simulator.do_circuit(sub_circuit)
+            # correct msmt values if ancilla feedback
+            if ancilla_fb: 
+                for index, qubit in enumerate([q.value for q in  ts]):
+                    msmt = simulator.current_measurement_record()[index-8]
+                    if qubit in z_ancillas:
+                        if z_ancillas_dict[qubit]['msmt_last'] == None:  
+                            z_ancillas_dict[qubit]['msmt_last'] = [msmt]
+                        else: 
+                            z_ancillas_dict[qubit]['msmt_last'].append(msmt)
+                            # if len(z_ancillas_dict[qubit]['msmt_last']) > 2:
+                            if len(z_ancillas_dict[qubit]['msmt_last']) > 3:
+                                z_ancillas_dict[qubit]['msmt_last'].pop(0) 
+                    elif qubit in x_ancillas: 
+                        if x_ancillas_dict[qubit]['msmt_last'] == None:  
+                            x_ancillas_dict[qubit]['msmt_last'] = [msmt]
+                        else: 
+                            x_ancillas_dict[qubit]['msmt_last'].append(msmt)
+                            # if len(x_ancillas_dict[qubit]['msmt_last']) > 2:
+                            if len(x_ancillas_dict[qubit]['msmt_last']) > 3:
+                                x_ancillas_dict[qubit]['msmt_last'].pop(0) 
+        
+        elif instruction.name == 'DEPOLARIZE2' :
+            ts = instruction.targets_copy()
+            qubits = [q.value for q in ts]
+            for j in range(int(len(qubits)/2)):
+                simulator.depolarize2(qubits[2*j], qubits[2*j + 1], p = after_clifford_depolarization_2)
+            simulator.do_circuit(sub_circuit)
+        else: 
+            simulator.do_circuit(sub_circuit)
+        
+    return simulator
+        
+import pymatching
+
+def count_logical_errors(circuit: stim.Circuit,   num_shots: int,
+                         p_amp_ge = 0.01, ancilla_fb = False, after_clifford_depolarization_2 =0.01) -> int:
+    # Sample the circuit.
+    # sampler = circuit.compile_detector_sampler()
+    # detection_events, observable_flips = sampler.sample(num_shots, separate_observables=True)
+    
+    
+    # surface_code_circuit.detector_error_model(decompose_errors=True)
+    measurements = []
+    for i in range(num_shots):
+        simulator = simulate_circuit(circuit,p_amp_ge=p_amp_ge, ancilla_fb = ancilla_fb, \
+            after_clifford_depolarization_2 =after_clifford_depolarization_2 )
+        measurements.append(simulator.current_measurement_record())
+    measurements = np.array(measurements, dtype=np.bool_)
+    # print('measurements',measurements)
+    detection_events, observable_flips = circuit.compile_m2d_converter().convert(measurements=measurements, separate_observables = True, append_observables=False) #, 
+
+    # Configure a decoder using the circuit.
+    detector_error_model = circuit.detector_error_model(decompose_errors=True)
+    matcher = pymatching.Matching.from_detector_error_model(detector_error_model)
+
+    # Run the decoder.
+    predictions = matcher.decode_batch(detection_events)
+
+    # Count the mistakes.
+    num_errors = 0
+    for shot in range(num_shots):
+        actual_for_shot = observable_flips[shot]
+        predicted_for_shot = predictions[shot]
+        if not np.array_equal(actual_for_shot, predicted_for_shot):
+            num_errors += 1
+    return num_errors
+
+import scipy
+def plot_data_S17(decoding_rounds = 3,  
+                  after_clifford_depolarization= 0.01, 
+                  after_clifford_depolarization_2 = 0.01,
+                  after_reset_flip_probability = 0.00, 
+                  before_measure_flip_probability = 0.00,  
+                  before_round_data_depolarization=0.01, 
+                  n_points = 50, max_err = 0.5, runs = 1000, plot = True, 
+                  ancilla_fb = False):
+    surface_code_circuit = stim.Circuit.generated(
+        "surface_code:rotated_memory_z",
+        rounds=decoding_rounds,
+        distance=3,
+        after_clifford_depolarization= after_clifford_depolarization,
+        after_reset_flip_probability=after_reset_flip_probability,
+        before_measure_flip_probability= before_measure_flip_probability,
+        before_round_data_depolarization=before_round_data_depolarization
+    )
+
+    if n_points ==1: 
+        p_amp_ge_array = np.linspace(max_err, max_err, 1)  # only simulate single point. 
+    else: 
+        p_amp_ge_array = np.linspace(0.00, max_err, n_points )
+    logical_err_array = []
+
+    for p_amp_ge in p_amp_ge_array: 
+        n_errors = count_logical_errors(surface_code_circuit,  
+                                        runs, 
+                                        p_amp_ge=p_amp_ge, 
+                                        ancilla_fb = ancilla_fb, 
+                                        after_clifford_depolarization_2 = after_clifford_depolarization_2 )
+        logical_err = n_errors/runs
+        logical_err_array.append(logical_err)
+        
+    if plot:    
+        param, vars = scipy.optimize.curve_fit(lambda t,a,b: a*np.exp(b*t),  p_amp_ge_array,  logical_err_array)
+        if ancilla_fb:
+            plt.plot(p_amp_ge_array,logical_err_array, label='raw: fb', color='r', marker='o')
+            plt.plot(p_amp_ge_array, param[0]*np.exp(param[1]*p_amp_ge_array),label= f'fit: fb (A:{ param[0]}, B: {param[1]})',color='r', marker = 'x' )
+        else: 
+            plt.plot(p_amp_ge_array,logical_err_array, label='raw: no fb', color='b', marker='o')
+            plt.plot(p_amp_ge_array, param[0]*np.exp(param[1]*p_amp_ge_array),label= f'fit: no fb  (A:{ param[0]}, B: {param[1]})',color='b', marker = 'x' )
+            
+        plt.legend()
+        plt.title(f' S17 Logical Z-error on  vs. p_eg on ancillas (dec. rnds: {decoding_rounds}, depol gate err: {after_clifford_depolarization}, data_err before rnd: {before_round_data_depolarization}, runs: {runs})')
+        plt.xlabel('p_eg ')
+        plt.ylabel('p_error')
+
+    return p_amp_ge_array, logical_err_array
+
+
+
+
+def plot_data_S17_rounds(decoding_rounds = 3,  
+                  after_clifford_depolarization= 0.001, 
+                  after_clifford_depolarization_2 = 0.01,
+                  after_reset_flip_probability = 0.001, 
+                  before_measure_flip_probability = 0.001,  
+                  before_round_data_depolarization=0.01, 
+                  n_points = 50, max_err = 0.5, runs = 1000, plot = True, 
+                  ancilla_fb = False):
+    
+
+    logical_err_array = []
+    
+    rounds_array = range(3, decoding_rounds, 3)
+    for rounds in rounds_array:
+        surface_code_circuit = stim.Circuit.generated(
+            "surface_code:rotated_memory_z",
+            rounds=rounds,
+            distance=3,
+            after_clifford_depolarization= after_clifford_depolarization,
+            after_reset_flip_probability=after_reset_flip_probability,
+            before_measure_flip_probability= before_measure_flip_probability,
+            before_round_data_depolarization=before_round_data_depolarization
+        )
+
+        n_errors = count_logical_errors(surface_code_circuit,  
+                                        runs, 
+                                        p_amp_ge=max_err, 
+                                        ancilla_fb = ancilla_fb, 
+                                        after_clifford_depolarization_2= after_clifford_depolarization_2)
+        logical_err = n_errors/runs
+        logical_err_array.append(logical_err)
+        
+    if plot:    
+        param, vars = scipy.optimize.curve_fit(lambda t,b: 1-np.exp(-b*t),  rounds_array,  logical_err_array)
+        print(param)
+        if max_err == 0 and not ancilla_fb:
+            plt.plot(rounds_array,logical_err_array, label= rf'raw: $p_{{ge}}:${max_err}', color='r', marker='o')
+            plt.plot(rounds_array, 1-np.exp(-param[0]*rounds_array),label= rf'fit: $ p_{{ge}}:${max_err} ($T_1$:{ round(1/param[0],1)})',\
+                linestyle = 'dotted', color='r', marker = 'x' )
+
+        elif max_err == 0 and ancilla_fb:
+            plt.plot(rounds_array,logical_err_array, label=rf'raw: fb, $p_{{ge}}:${max_err}', color='g', marker='o')
+            plt.plot(rounds_array,1-np.exp(-param[0]*rounds_array),label= rf'fit: fb, $p_{{ge}}:${max_err} ($T_1$:{ round(1/param[0],1)})' ,\
+                linestyle = 'dotted',color='g',marker = 'x' )
+             
+        elif ancilla_fb:
+            plt.plot(rounds_array,logical_err_array, label=rf'raw: fb, $p_{{ge}}:${max_err}', color='b', marker='o')
+            plt.plot(rounds_array, 1-np.exp(-param[0]*rounds_array),label= rf'fit: fb, $p_{{ge}}:${max_err} ($T_1$:{ round(1/param[0],1)})',\
+                linestyle = 'dotted',color='b',marker = 'x' )
+            
+        else: 
+            plt.plot(rounds_array,logical_err_array, label=rf'raw: $p_{{ge}}:${max_err}', color='k', marker='.')
+            plt.plot(rounds_array, 1-np.exp(-param[0]*rounds_array),label= rf'fit: $p_{{ge}}:${max_err} ($T_1$:{ round(1/param[0],1)})',\
+                linestyle = 'dotted',color='k',marker = 'x' )
+            
+        plt.legend(  loc='upper left', borderaxespad=0)
+        # plt.title(f' S17 Logical Z-error on  vs. p_eg on ancillas (dec. rnds: {decoding_rounds}, depol gate err: {after_clifford_depolarization}, p_eg: {max_err} \
+        #     ,data_err before rnd: {before_round_data_depolarization}, after_reset_flip_probability; {after_reset_flip_probability},before_measure_flip_probability:\
+        #         {before_measure_flip_probability} runs: {runs})')
+        plt.title(rf' S17, logical error probability $E_L$ (for $|0\rangle_L$), $p_{{ge}}$: {max_err}')
+        plt.xlabel('rounds')
+        plt.ylabel('logical error probability $E_L$ ')
+
+    return rounds_array, logical_err_array
\ No newline at end of file