Qiskit statevector simulator returns wrong result when skipping classical bits.

[nathanieltornow]

Hi!

I noticed that when "skipping" a classical bit, meaning not measuring to some bit $x$ but to bit $y$ with $ind(x) < ind(y)$, the cusvaer-enabled statevector simulator returns the wrong result, simply ignoring the result on bit $y$. This does not happen when not enabling cusvaer.

See the example below:

from qiskit.circuit import QuantumCircuit
from qiskit_aer import Aer

sim = Aer.get_backend("aer_simulator_statevector")
sim.set_option("cusvaer_enable", True)

circuit = QuantumCircuit(4, 5)
circuit.h(0)
circuit.cx(range(3), range(1, 4))
circuit.measure([0, 1, 2, 3], [0, 1, 2, 4])

print(circuit)

wrong_result = sim.run(circuit, shots=10000).result().get_counts()
print("Wrong:", wrong_result)

sim.set_option("cusvaer_enable", False)

correct_result = sim.run(circuit, shots=10000).result().get_counts()
print("Correct:", correct_result)

     ┌───┐          ┌─┐              
q_0: ┤ H ├──■───────┤M├──────────────
     └───┘┌─┴─┐     └╥┘     ┌─┐      
q_1: ─────┤ X ├──■───╫──────┤M├──────
          └───┘┌─┴─┐ ║      └╥┘┌─┐   
q_2: ──────────┤ X ├─╫───■───╫─┤M├───
               └───┘ ║ ┌─┴─┐ ║ └╥┘┌─┐
q_3: ────────────────╫─┤ X ├─╫──╫─┤M├
                     ║ └───┘ ║  ║ └╥┘
c: 5/════════════════╩═══════╩══╩══╩═
                     0       1  2  4 
Wrong: {'00000': 4993, '00111': 5007}
Correct: {'10111': 5026, '00000': 4974}

I suspect there might be some missing bookkeeping of the indices of classical registers in the cusvaer backend.

Thanks!
Nate

[yangcal]

Thanks for reporting. @ymagchi can you take a look?

[nathanieltornow]

Perhaps, something like this could be done on the frontend:

def clbit_mapping(circuit: QuantumCircuit) -> tuple[QuantumCircuit, list[int]]:
    active_clbit_indices = sorted(
        set(circuit.clbits.index(clbit) for instr in circuit for clbit in instr.clbits),
    )
    mapping = {clbit: i for i, clbit in enumerate(active_clbit_indices)}

    new_creg = ClassicalRegister(len(active_clbit_indices), name="c")
    new_circuit = QuantumCircuit(*circuit.qregs, new_creg)
    for instr in circuit:
        new_circuit.append(
            instr.operation,
            instr.qubits,
            [mapping[circuit.clbits.index(clbit)] for clbit in instr.clbits],
        )

    return new_circuit, active_clbit_indices


def remap_results(bitstr: str, num_clbits: int, active_indices: list[int]) -> str:
    assert len(bitstr) <= num_clbits
    bits = ["0"] * num_clbits
    for i, idx in enumerate(active_indices):
        bits[-idx - 1] = bitstr[i]
    return "".join(bits)


circuit2, active_indices = clbit_mapping(circuit)
corrected_res = sim.run(circuit2, shots=10000).result().get_counts()
print("corrected:", corrected_res)
print({remap_results(k, circuit.num_clbits, active_indices): v for k, v in corrected_res.items()})

[ymagchi]

@nathanieltornow Thank you so much for reporting this. I confirmed that it is reproducible with 24.03 images. We will work on its fix.