# Qiskit DensityMatrix.from_instruction when snapshots are present

I have a quantum circuit in which I apply snapshots like this during setup:

    for qubit_index in range(0, num_qubits):

if qubit_index < num_qubits-1:
next_qubit_index = qubit_index + 1
snapshot_name_string = "snapshot_" + str(qubit_index) + "_" + str(next_qubit_index)
print("--- Setting " + snapshot_name_string)
quantum_circuit.snapshot(snapshot_name_string, qubits = [qubit_index,next_qubit_index])



This so I can compute the amplitudes of sub-systems. For example, if I have 10 qubits, but want to get probabilities of only the 3rd and fourth qubit being in their states (00, 01, 10, 11), I use:

result = execute(quantum_circuit, backend, shots=1).result()
snapshots = result.data()['snapshots']['statevector'][snapshot_name_string][0].tolist()


I need the probabilities to correctly sample projective measurements before applying Clifford gates to (qubit_index,next_qubit_index) during the evolution. However, later on in the evolution when I try to analyze the density matrix:

rho = qi.DensityMatrix.from_instruction(quantum_circuit)
# just note I then go on to compute the partial trace
reduced_rho = qi.partial_trace(rho, subsystem_range)


I now get an error QiskitError: 'Cannot apply Instruction: snapshot' from qi.DensityMatrix.from_instruction(). The density matrix from_instruction was working fine before I implemented the snapshots (when I was randomly selected a projective measurement based on uniform distribution rather than based on amplitudes).

Does anyone know how I can reconcile the DensityMatrix with the snapshots? I don't see a DensityMatrix parameter that just says 'use the whole system instead of a snapshot.' Also if their is an easy/better way to get amplitudes of sub-system states without snapshots that would be perfectly acceptable answer as well!

<thank|you!>

Since you already executed the circuit and get the snapshots you want, a workaround is to remove the snapshot instructions before calling DensityMatrix.from_instruction:

from itertools import count
for index, instruction in zip(count(len(circ.data) - 1, -1), reversed(circ.data)):
print(index, instruction[0].name)
if instruction[0].name == 'snapshot':
del circ.data[index]

DensityMatrix.from_instruction(circ)


You may also replace them with Barrier instructions or I gates:

from qiskit.circuit import Barrier
quantum_circuit.data = [(Barrier(_inst[0].num_qubits), _inst[1], _inst[2]) if _inst[0].name == 'snapshot' else _inst for _inst in quantum_circuit.data]

DensityMatrix.from_instruction(circ)