Syndrome showing unexpected results for the [5,1,3] QEC code in qiskit

Question

I have used the encoding procedure mentioned on page 35 of https://arxiv.org/pdf/quant-ph/9705052.pdf for the [5,1,3] code and then built the circuit for measuring the stabilizers (XZZXI, IXZZX, XIXZZ, ZXIXZ). I have used the respective Hadamard transformations to convert the X measurements to the Z measurements and by connecting it with the ancilla qubits through the respective gates I constructed the circuit to measure the syndrome. But instead of showing '0000' as the syndrome(since I have included no errors), it is showing results like,

{'1111 00000': 628, '1011 00000': 659, '0100 00000': 616, '1001 00000': 626, '0000 00000': 648, '1101 00000': 653, '0001 00000': 657, '1010 00000': 641, '0011 00000': 598, '0101 00000': 624, '0010 00000': 649, '0110 00000': 601, '1100 00000': 643, '0111 00000': 694, '1110 00000': 637, '1000 00000': 666}

where the first four qubits are the ancilla qubits. To make matters simple I have used the encoded state directly(so that there is no fault in encoding). Could there be a problem in the syndrome measurement circuit or it has something to do with the way qiskit's qasm_simulator?

I have attached my code below,

from qiskit import *
import numpy as np
from sympy.physics.quantum.tensorproduct import  TensorProduct

es=list(range(0,32))

import math
a=1/math.sqrt(2)
es[0]=es[18]=es[9]=es[20]=es[10]=a/4
es[27]=es[6]=es[24]=es[29]=es[3]=es[30]=es[15]=es[17]=es[12]=es[23]=-a/4
es[5]=a/4

    

b=1/math.sqrt(2)
es[31]=es[13]=es[22]=es[11]=es[21]=b/4
es[4]=es[25]=es[7]=es[2]=es[28]=es[1]=es[16]=es[14]=es[19]=es[8]=-b/4
es[26]=b/4

ph_qr=QuantumRegister(5,"physical_qubits")
ph_cr=ClassicalRegister(5,"measured_physical_bits")
an_qr=QuantumRegister(4,"ancilla_qubits")
an_cr=ClassicalRegister(4,"measured_ancilla_bits")
syndrome=QuantumCircuit(ph_qr, an_qr, ph_cr, an_cr)
syndrome.initialize(es, ph_qr)

syndrome_nf=syndrome.copy()
syndrome_nf.h(ph_qr[2])
syndrome_nf.h(ph_qr[3])
syndrome_nf.h(an_qr[0])
syndrome_nf.h(an_qr[1])
syndrome_nf.h(an_qr[2])
syndrome_nf.h(an_qr[3])
syndrome_nf.cz(ph_qr[4], an_qr[0])
syndrome_nf.cz(ph_qr[3], an_qr[0])
syndrome_nf.h(ph_qr[3])
syndrome_nf.cz(ph_qr[2], an_qr[0])
syndrome_nf.h(ph_qr[2])
syndrome_nf.cz(ph_qr[1], an_qr[0])
syndrome_nf.h(ph_qr[4])
syndrome_nf.cz(ph_qr[4], an_qr[1])
syndrome_nf.h(ph_qr[4])
syndrome_nf.h(ph_qr[3])
syndrome_nf.cz(ph_qr[3], an_qr[1])
syndrome_nf.h(ph_qr[3])
syndrome_nf.cz(ph_qr[2], an_qr[1])
syndrome_nf.cz(ph_qr[0], an_qr[1])
syndrome_nf.h(ph_qr[4])
syndrome_nf.cz(ph_qr[4], an_qr[2])
syndrome_nf.h(ph_qr[4])
syndrome_nf.cz(ph_qr[3], an_qr[2])
syndrome_nf.cz(ph_qr[1], an_qr[2])
syndrome_nf.h(ph_qr[0])
syndrome_nf.cz(ph_qr[0] ,an_qr[2])
syndrome_nf.cz(ph_qr[4], an_qr[3])
syndrome_nf.cz(ph_qr[2], an_qr[3])
syndrome_nf.h(ph_qr[1])
syndrome_nf.cz(ph_qr[1], an_qr[3])
syndrome_nf.h(ph_qr[1])

syndrome_nf.cz(ph_qr[0], an_qr[3])
syndrome_nf.h(ph_qr[0])
syndrome_nf.h(an_qr[0])
syndrome_nf.h(an_qr[1])
syndrome_nf.h(an_qr[2])
syndrome_nf.h(an_qr[3])
syndrome_nf.draw()

syndrome_nf.measure(an_qr, an_cr)

simulator=qiskit.Aer.get_backend('qasm_simulator')
result=execute(syndrome_nf, backend=simulator, shots=10240).result()
count=result.get_counts()

print(count)

Answer 1

Look at the start of the multiqubit section in this tutorial. In particular the section on basis vector ordering. I found the ordering of qubits to be very strange in qiskit possibly this is your error as well? For example the state $|10\rangle$ corresponds to qubit 0 being in state $|0\rangle$ and qubit 1 in state $|1\rangle$, contrary to what you might expect.