Add psi exchange, start writing README.md

Author: Randl

README.md

@@ -6,7 +6,26 @@ in Qiskit.
 
 ## GS preparation
 ### Matching boundary conditions
+For matching boundary condition, boundary plaquettes are all of the same type and the ground state is unique.
+We implement linear algorithm for preparing ground state of the toric code proposed in the paper.
 ### Mixed boundary conditions
+For mixed boundary condition, boundary plaquettes are of different types, and there is ground state degeneracy,
+which allows us to encode logical qubits in the system state.
+Not implemented yet.
+
 ## Entropy and topological entropy
+We implement the measurement of the second Rényi entropy as described in paper, and use it
+to calculate topological entropy, acquiring non-trivial value for various subsystems.
+
+For 2x2 and 2x3 subsystem it is possible to perform determenistic calculation, while for 
+3x3 system only randomized calculation is feasible.
+
+## Braiding
+There are 4 particle types in toric code: `e`, `m`, `psi` and `1`, for total of 6 possible 
+mutual statistics and 3 exchange statistic. Out of those 4 are non-trivial -- `em`, `epsi`,
+`mpsi` and `psipsi`, resulting in phase of π.
 
-## Braiding
+We implement the required operators (without optiimization), and demonstrate part of the braidings
+and exchanges.
+## Logical qubit.
+Not implemented yet.

tests/test_braiding.py

@@ -4,7 +4,7 @@
 from qiskit import Aer
 from qiskit import transpile
 
-from topo_braiding import create_e_particles, create_m_particles, apply_cxxxx_on_square
+from topo_braiding import create_e_particles, create_m_particles, apply_cxxxx_on_square, apply_cxxyyzz_on_rectangle
 from toric_code import get_toric_code
 
 
@@ -41,6 +41,76 @@ def test_em_braiding(self):
             expected_res = -1. if sq in [(0, 2), (0, 3)] else 1.
             np.testing.assert_allclose(cos_theta, expected_res)
 
+    def test_psipsi_exchange(self):
+        x, y = 5, 7
+        backend = Aer.get_backend('aer_simulator')
+        sq = (0, 3)
+        tc = get_toric_code(x, y, classical_bit_count=1, ancillas_count=1)
+
+        x_strings = [[(2, 1), (2, 2)],
+                     [(1, 4), (3, 4)]]
+        for x_string in x_strings:
+            create_e_particles(tc, x_string)
+
+        z_strings = [[(2, 3), (3, 3), (4, 3)],
+                     [(3, 1), (4, 1), (3, 2)], ]
+        for z_string in z_strings:
+            create_m_particles(tc, z_string)
+
+        tc.circ.h(tc.ancillas[0])
+        apply_cxxyyzz_on_rectangle(tc, sq)
+        tc.circ.h(tc.ancillas[0])
+
+        tc.circ.measure(tc.ancillas[0], 0)
+        Nshots = 10000
+        job = backend.run(transpile(tc.circ, backend), shots=Nshots)
+        result = job.result()
+        counts = result.get_counts(tc.circ)
+
+        if '0' not in counts:
+            counts['0'] = 0
+        if '1' not in counts:
+            counts['1'] = 0
+        print(counts)
+
+        cos_theta = (counts['0'] - counts['1']) / Nshots
+        expected_res = -1.
+        np.testing.assert_allclose(cos_theta, expected_res)
+
+    def test_psione_braiding(self):
+        x, y = 5, 7
+        backend = Aer.get_backend('aer_simulator')
+        sq = (0, 3)
+        tc = get_toric_code(x, y, classical_bit_count=1, ancillas_count=1)
+
+        x_strings = [[(1, 4), (3, 4)]]
+        for x_string in x_strings:
+            create_e_particles(tc, x_string)
+
+        z_strings = [[(2, 3), (3, 3), (4, 3)], ]
+        for z_string in z_strings:
+            create_m_particles(tc, z_string)
+
+        tc.circ.h(tc.ancillas[0])
+        apply_cxxyyzz_on_rectangle(tc, sq)
+        tc.circ.h(tc.ancillas[0])
+
+        tc.circ.measure(tc.ancillas[0], 0)
+        Nshots = 10000
+        job = backend.run(transpile(tc.circ, backend), shots=Nshots)
+        result = job.result()
+        counts = result.get_counts(tc.circ)
+
+        if '0' not in counts:
+            counts['0'] = 0
+        if '1' not in counts:
+            counts['1'] = 0
+        print(counts)
+
+        cos_theta = (counts['0'] - counts['1']) / Nshots
+        expected_res = 1.
+        np.testing.assert_allclose(cos_theta, expected_res)
+
     def test_e_trivial_braiding(self):
         x, y = 5, 7
         backend = Aer.get_backend('aer_simulator')

tests/test_entropy.py

@@ -6,7 +6,7 @@
 from qiskit import transpile
 from tqdm import tqdm
 
-from topo_entropy import ABC_DIVISION_2x2, ABC_DIVISION_2x3_LEFT, ABC_DIVISION_2x3_RIGHT,ABC_DIVISION_3x3
+from topo_entropy import ABC_DIVISION_2x2, ABC_DIVISION_2x3_LEFT, ABC_DIVISION_2x3_RIGHT, ABC_DIVISION_3x3
 from topo_entropy import calculate_s_subsystems
 from topo_entropy import get_all_2x2_non_corner, get_all_2x3_non_corner, get_all_3x3_non_corner
 from topo_entropy import get_all_2x3_left_non_corner, get_all_2x3_right_non_corner
@@ -30,7 +30,7 @@ def test_topo_entropy(backend, size, qubits, subsystems, expected_values, type='
         elif type == 'pauli':
             tc.measure_pauli(qubits, gates)
 
-        job = backend.run(transpile(tc.circ, backend), shots=15000) # note: number of shots is important
+        job = backend.run(transpile(tc.circ, backend), shots=15000)  # note: number of shots is important
         result = job.result()
         counts = result.get_counts(tc.circ)
         all_counts.append(counts)
@@ -71,27 +71,30 @@ def test_2x3_entropy_pauli(self):
         backend_sim = Aer.get_backend('aer_simulator')
         x, y = 5, 7
         for qubits in get_all_2x3_left_non_corner((x, y)):
-            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_LEFT, expected_values, type='pauli')
+            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_LEFT, expected_values, type='pauli',
+                              rtol=0.03)
         for qubits in get_all_2x3_right_non_corner((x, y)):
-            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_RIGHT, expected_values, type='pauli')
+            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_RIGHT, expected_values, type='pauli',
+                              rtol=0.03)
 
     def test_2x3_entropy_haar(self):
         expected_values = [(2., 2., 2.), (4., 4., 3.), (4.,)]
 
         backend_sim = Aer.get_backend('aer_simulator')
         x, y = 5, 7
+        for qubits in get_all_2x3_right_non_corner((x, y)):
+            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_RIGHT, expected_values, type='haar',
+                              cnt=100)
         for qubits in get_all_2x3_left_non_corner((x, y)):
             test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_LEFT, expected_values, type='haar', cnt=100)
-        for qubits in get_all_2x3_right_non_corner((x, y)):
-            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_2x3_RIGHT, expected_values, type='haar', cnt=100)
 
     def test_3x3_entropy_haar(self):
         expected_values = [(3., 3., 3.), (6., 5., 4.), (5.,)]
 
         backend_sim = Aer.get_backend('aer_simulator')
         x, y = 5, 7
         for qubits in get_all_3x3_non_corner((x, y)):
-            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_3x3, expected_values,  type='haar', cnt=1000)
+            test_topo_entropy(backend_sim, (x, y), qubits, ABC_DIVISION_3x3, expected_values, type='haar', cnt=1000)
 
 
 if __name__ == '__main__':

topo_braiding.py

@@ -1,8 +1,8 @@
 import numpy as np
 from qiskit import transpile
-from qiskit.circuit.library import MCXGate
 
 from toric_code import get_toric_code
+from toric_code_matching import is_inside_matching
 
 sz = np.array([[1, 0], [0, -1]])
 sx = np.array([[0, 1], [1, 0]])
@@ -26,21 +26,46 @@ def apply_cxxxx_on_square(tc, upper_corner):
         locs = [(x, y), (x + 1, y), (x + 2, y), (x + 1, y + 1)]
     else:
         locs = [(x, y), (x + 1, y - 1), (x + 2, y), (x + 1, y)]
-    tc.circ.mct(control_qubits=[tc.ancillas[0]], target_qubit=[tc.regs[l[0]][l[1]] for l in locs])
 
+    if not all([is_inside_matching((tc.x, tc.y), l) for l in locs]):
+        return
+    tc.circ.mct(control_qubits=[tc.ancillas[0]], target_qubit=[tc.regs[l[0]][l[1]] for l in locs])
 
-def apply_cxxyyzz_on_square(tc, upper_corner):
-    # cXXXX = np.kron(np.array([[1, 0], [0, 0]]), np.kron(sx, np.kron(sx, np.kron(sx, sx)))) + \
-    #         np.kron(np.array([[0, 0], [0, 1]]), np.kron(s0, np.kron(s0, np.kron(s0, s0))))
 
+def apply_czzz_on_square(tc, upper_corner):
+    cZZZZ = np.kron(np.array([[1, 0], [0, 0]]), np.kron(sz, np.kron(sz, np.kron(sz, sz)))) + \
+            np.kron(np.array([[0, 0], [0, 1]]), np.kron(s0, np.kron(s0, np.kron(s0, s0))))
     x, y = upper_corner
     if x % 2 == 0:
         locs = [(x, y), (x + 1, y), (x + 2, y), (x + 1, y + 1)]
     else:
         locs = [(x, y), (x + 1, y - 1), (x + 2, y), (x + 1, y)]
 
-    tc.circ.mct(control_qubits=[tc.ancillas[0]], target_qubit=[tc.regs[l[0]][l[1]] for l in locs])
-    # tc.circ.unitary(cXXXX, [tc.regs[l[0]][l[1]] for l in locs] + [tc.ancillas[0]])
+    if not all([is_inside_matching((tc.x, tc.y), l) for l in locs]):
+        return
+
+    tc.circ.unitary(cZZZZ, [tc.regs[l[0]][l[1]] for l in locs] + [tc.ancillas[0]])
+
+
+def apply_cxxyyzz_on_rectangle(tc, upper_corner, left=True):
+    cXXYYZZ = np.kron(np.array([[1, 0], [0, 0]]), np.kron(sx, np.kron(sx, np.kron(sy, np.kron(sy, np.kron(sz, sz)))))) + \
+              np.kron(np.array([[0, 0], [0, 1]]), np.kron(s0, np.kron(s0, np.kron(s0, np.kron(s0, np.kron(s0, s0))))))
+
+    x, y = upper_corner
+    if left:
+        if x % 2 == 0:
+            locs = [(x, y), (x + 1, y + 1), (x + 1, y), (x + 2, y), (x + 2, y - 1), (x + 3, y)]
+        else:
+            locs = [(x, y), (x + 1, y), (x + 1, y - 1), (x + 2, y), (x + 1, y - 1), (x + 2, y - 1)]
+    else:
+        if x % 2 == 0:
+            locs = [(x, y), (x + 1, y), (x + 1, y + 1), (x + 2, y), (x + 2, y + 1), (x + 3, y + 1)]
+        else:
+            locs = [(x, y), (x + 1, y - 1), (x + 1, y), (x + 2, y), (x + 2, y + 1), (x + 3, y)]
+
+    if not all([is_inside_matching((tc.x, tc.y), l) for l in locs]):
+        return
+    tc.circ.unitary(cXXYYZZ, [tc.regs[l[0]][l[1]] for l in locs[::-1]] + [tc.ancillas[0]])
 
 
 def em_braiding_phase(backend, x, y):

topo_entropy.py

@@ -10,7 +10,7 @@
 
 ABC_DIVISION_2x2 = [(0, 1), (2,), (3,)]
 ABC_DIVISION_2x3_RIGHT = [(1, 3), (0, 2), (4, 5)]
-ABC_DIVISION_2x3_LEFT = [(2, 4), (0, 1),  (3, 5)]
+ABC_DIVISION_2x3_LEFT = [(2, 4), (0, 1), (3, 5)]
 ABC_DIVISION_3x3 = [(0, 1, 3), (2, 5, 7), (4, 6, 8)]
 
 
@@ -122,6 +122,8 @@ def get_all_2x3_left_non_corner(size):
                     [is_corner_matching(size, s) for s in sys_l]):
                 all_sys.append(sys_l)
     return all_sys
+
+
 def get_all_2x3_right_non_corner(size):
     x, y = size
     all_sys = []
@@ -135,6 +137,8 @@ def get_all_2x3_right_non_corner(size):
                     [is_corner_matching(size, s) for s in sys_r]):
                 all_sys.append(sys_r)
     return all_sys
+
+
 def get_all_2x3_non_corner(size):
     return get_all_2x3_left_non_corner(size) + get_all_2x3_right_non_corner(size)
 
@@ -147,8 +151,7 @@ def get_all_3x3_non_corner(size):
             if rx % 2 == 0:
                 sys = [(-1, -1)]  # skip, why
             else:
-                sys = (rx, ry), (rx + 1, ry - 1), (rx + 1, ry), (rx + 2, ry - 1), (rx + 2, ry), (rx + 2, ry + 1), (
-                rx + 3, ry - 1), (rx + 3, ry), (rx + 4, ry),
+                sys = (rx, ry), (rx + 1, ry - 1), (rx + 1, ry), (rx + 2, ry - 1), (rx + 2, ry), (rx + 2, ry + 1), (rx + 3, ry - 1), (rx + 3, ry), (rx + 4, ry),
             if all([is_inside_matching(size, s) for s in sys]) and not any([is_corner_matching(size, s) for s in sys]):
                 all_sys.append(sys)
     return all_sys

toric_code.py

@@ -5,6 +5,6 @@
 def get_toric_code(x, y, classical_bit_count=4, ancillas_count=0, boundary_condition='matching'):
     assert boundary_condition in ('mixed', 'matching')
     if boundary_condition == 'matching':
-        return ToricCodeMatching(x, y, classical_bit_count,ancillas_count)
+        return ToricCodeMatching(x, y, classical_bit_count, ancillas_count)
     elif boundary_condition == 'mixed':
-        return ToricCodeMixed(x, y, classical_bit_count,ancillas_count)
+        return ToricCodeMixed(x, y, classical_bit_count, ancillas_count)

toric_code_matching.py

@@ -83,10 +83,9 @@ def __init__(self, x, y, classical_bit_count=4, ancillas_count=0):
         self.star_x, self.star_y = self.x, self.y // 2 + 1
         # print(self.plaquette_x, self.plaquette_y)
         self.regs = [QuantumRegister(self.x - 1, f'l{lev}') if lev % 2 == 0 else QuantumRegister(self.x, f'l{lev}')
-                     for
-                     lev in range(self.y)]  # first coordinate is row index, second is column index
+                     for lev in range(self.y)]  # first coordinate is row index, second is column index
         self.c_reg = ClassicalRegister(classical_bit_count)
-        if ancillas_count>0:
+        if ancillas_count > 0:
             self.ancillas = QuantumRegister(ancillas_count)
             self.circ = QuantumCircuit(*self.regs, self.ancillas, self.c_reg)
         else: