Fix snapshot density matrix (#374)

* change density matrix snapshot from singleshot to average snapshot

* adds State name to error message for invalid circuit instructions

* add density matrix snapshot tests

CHANGELOG.md

@@ -23,9 +23,11 @@ Added
 - Added tests for the Fredkin gate (#357)
 - Added tests for the cu1 gate (#360)
 - Added tests for statevector and stabilizer snapshots (\#355)
+- Added tests for density matrix snapshot (\#374)
 
 Changed
 -------
+- Change signature of density matrix snapshot extension (\#374)
 - Stabilizer snapshot returns stabilizer instead of full Clifford table (\#355)
 - Signature of SnapshotStatevector and SnapshotStabilizer (\#355)
 - Changed all names from tensor_network_state to matrix_product_state (\#356)

qiskit/providers/aer/extensions/snapshot_density_matrix.py

@@ -21,38 +21,41 @@
 class SnapshotDensityMatrix(Snapshot):
     """Snapshot instruction for density matrix method of Qasm simulator."""
 
-    def __init__(self,
-                 label,
-                 num_qubits=0,
-                 num_clbits=0,
-                 params=None):
+    def __init__(self, label, num_qubits):
+        """Create a density matrix state snapshot instruction.
 
-        super().__init__(label, 'density_matrix', num_qubits, num_clbits, params)
+        Args:
+            label (str): the snapshot label.
+            num_qubits (int): the number of qubits to snapshot.
 
+        Raises:
+            ExtensionError: if snapshot is invalid.
+        """
+        super().__init__(label,
+                         snapshot_type='density_matrix',
+                         num_qubits=num_qubits)
+
+
+def snapshot_density_matrix(self, label, qubits=None):
+    """Take a density matrix snapshot of simulator state.
 
-def snapshot_density_matrix(self,
-                            label,
-                            qubits=None,
-                            params=None):
-    """Take a density matrix snapshot of the internal simulator representation.
-    Works on all qubits, and prevents reordering (like barrier).
     Args:
         label (str): a snapshot label to report the result
-        qubits (list or None): the qubits to apply snapshot to [Default: None].
-        params (list or None): the parameters for snapshot_type [Default: None].
+        qubits (list or None): the qubits to apply snapshot to. If None all
+                               qubits will be snapshot [Default: None].
+
     Returns:
-        QuantumCircuit: with attached command
+        QuantumCircuit: with attached instruction.
+
     Raises:
-        ExtensionError: malformed command
+        ExtensionError: if snapshot is invalid.
     """
 
     snapshot_register = Snapshot.define_snapshot_register(self, label, qubits)
 
     return self.append(
-        SnapshotDensityMatrix(
-            label,
-            num_qubits=len(snapshot_register),
-            params=params), snapshot_register)
+        SnapshotDensityMatrix(label, num_qubits=len(snapshot_register)),
+        snapshot_register)
 
 
 QuantumCircuit.snapshot_density_matrix = snapshot_density_matrix

src/base/state.hpp

@@ -256,7 +256,8 @@ std::string State<state_t>::invalid_opset_message(const Operations::OpSet &opset
   // We can't print OpTypes so we add a note if there are invalid
   // instructions other than gates or snapshots
   if (bad_instr && (!bad_gates && !bad_snaps))
-    ss << " invalid non gate or snapshot instructions: opset={" << opset << "}";
+    ss << " invalid non gate or snapshot instructions in opset {" << opset << "}";
+  ss << " for " << name() << " method"; 
   return ss.str();
 }
 

src/simulators/densitymatrix/densitymatrix_state.hpp

@@ -443,7 +443,11 @@ void State<densmat_t>::apply_snapshot(const Operations::Op &op,
                                 op.name + "\'.");
   switch (it -> second) {
     case Snapshots::densitymatrix:
-      BaseState::snapshot_state(op, data, "density_matrix");
+      data.add_average_snapshot("density_matrix",
+                                op.string_params[0],
+                                BaseState::creg_.memory_hex(),
+                                BaseState::qreg_,
+                                false);
       break;
     case Snapshots::cmemory:
       BaseState::snapshot_creg_memory(op, data);
@@ -483,13 +487,16 @@ void State<densmat_t>::apply_snapshot(const Operations::Op &op,
 
 template <class densmat_t>
 void State<densmat_t>::snapshot_probabilities(const Operations::Op &op,
-                                               OutputData &data,
-                                               bool variance) {
+                                              OutputData &data,
+                                              bool variance) {
   // get probs as hexadecimal
   auto probs = Utils::vec2ket(measure_probs(op.qubits),
                               json_chop_threshold_, 16);
-  data.add_average_snapshot("probabilities", op.string_params[0],
-                            BaseState::creg_.memory_hex(), probs, variance);
+  data.add_average_snapshot("probabilities",
+                            op.string_params[0],
+                            BaseState::creg_.memory_hex(),
+                            probs,
+                            variance);
 }
 
 

src/simulators/qasm/qasm_controller.hpp

@@ -468,6 +468,18 @@ QasmController::simulation_method(const Circuit &circ,
                                   bool validate) const {
   // Check simulation method and validate state
   switch(simulation_method_) {
+    case Method::statevector: {
+      if (validate) {
+        if (simulation_precision_ == Precision::single_precision) {
+          Statevector::State<QV::QubitVector<float>> state; 
+          validate_state(state, circ, noise_model, true);
+        } else {
+          Statevector::State<QV::QubitVector<>> state; 
+          validate_state(state, circ, noise_model, true);
+        }
+      }
+      return Method::statevector;
+    }
     case Method::density_matrix: {
       if (validate)
         validate_state(DensityMatrix::State<>(), circ, noise_model, true);
@@ -504,7 +516,7 @@ QasmController::simulation_method(const Circuit &circ,
           validate_state(DensityMatrix::State<>(), circ, noise_model, false) &&
           check_measure_sampling_opt(circ, Method::density_matrix).first) {
         return Method::density_matrix;
-        }
+      }
       // Finally we check the statevector memory requirement for the
       // current number of qubits. If it fits in available memory we
       // default to the Statevector method. Otherwise we attempt to use
@@ -527,10 +539,16 @@ QasmController::simulation_method(const Circuit &circ,
     }
     // If we didn't select extended stabilizer above proceed to the default switch clause
     default: {
-      // Default method is statevector
+      // For default we use statevector followed by density matrix (for the case
+      // when the circuit contains invalid instructions for statevector)
+      if (validate_state(Statevector::State<>(), circ, noise_model, false)) {
+        return Method::statevector;
+      }
+      // If circuit contains invalid instructions for statevector throw a hail
+      // mary and try for density matrix.
       if (validate)
-        validate_state(Statevector::State<>(), circ, noise_model, true);
-      return Method::statevector;
+        validate_state(DensityMatrix::State<>(), circ, noise_model, true);
+      return Method::density_matrix;
     }
   }
 }

test/terra/backends/qasm_simulator/qasm_snapshot.py

@@ -30,6 +30,10 @@
     snapshot_state_pre_measure_statevector_nondeterministic,
     snapshot_state_post_measure_statevector_nondeterministic)
 
+def get_snapshots(data, label, snapshot_type):
+        """Format snapshots as list of Numpy arrays"""
+        return data.get("snapshots", {}).get(snapshot_type, {}).get(label, [])
+
 
 class QasmSnapshotStatevectorTests:
     """QasmSimulator snapshot statevector tests."""
@@ -40,14 +44,10 @@ class QasmSnapshotStatevectorTests:
     ]
     BACKEND_OPTS = {}
 
-    def statevector_snapshots(self, data, label):
+    @staticmethod
+    def statevector_snapshots(data, label):
         """Format snapshots as list of Numpy arrays"""
-        # Check snapshot entry exists in data
-        self.assertIn("snapshots", data)
-        self.assertIn("statevector", data["snapshots"])
-        self.assertIn(label, data["snapshots"]["statevector"])
-        # Format output as list of numpy array
-        snaps = data["snapshots"]["statevector"][label]
+        snaps = data.get("snapshots", {}).get("statevector", {}).get(label, [])
         statevecs = []
         for snap in snaps:
             svec = np.array(snap)
@@ -188,14 +188,10 @@ class QasmSnapshotStabilizerTests:
     SUPPORTED_QASM_METHODS = ['automatic', 'stabilizer']
     BACKEND_OPTS = {}
 
-    def stabilizer_snapshots(self, data, label):
-        """Format snapshots as list of Numpy arrays"""
-        # Check snapshot entry exists in data
-        self.assertIn("snapshots", data)
-        self.assertIn("stabilizer", data["snapshots"])
-        self.assertIn(label, data["snapshots"]["stabilizer"])
-        # Format output as list of numpy array
-        return data["snapshots"]["stabilizer"][label]
+    @staticmethod
+    def stabilizer_snapshots(data, label):
+        """Get stabilizer snapshots"""
+        return data.get("snapshots", {}).get("stabilizer", {}).get(label, [])
 
     @staticmethod
     def stabilizes_statevector(stabilizer, statevector):
@@ -340,3 +336,151 @@ def test_snapshot_stabilizer_post_measure_nondet(self):
                     stabilizer = snaps[j]
                     self.assertTrue(
                         self.stabilizes_statevector(stabilizer, statevec))
+
+
+class QasmSnapshotDensityMatrixTests:
+    """QasmSimulator snapshot density matrix tests."""
+
+    SIMULATOR = QasmSimulator()
+    SUPPORTED_QASM_METHODS = [
+        'automatic', 'density_matrix'
+    ]
+    BACKEND_OPTS = {}
+
+    @staticmethod
+    def density_snapshots(data, label):
+        """Format snapshots as list of Numpy arrays"""
+        # Check snapshot entry exists in data
+        snaps = data.get("snapshots", {}).get("density_matrix", {}).get(label, [])
+        # Convert nested lists to numpy arrays
+        output = {}
+        for snap_dict in snaps:
+            memory = snap_dict['memory']
+            mat = np.array(snap_dict['value'])
+            output[memory] = mat[:, :, 0] + 1j * mat[:, :, 1]
+        return output
+
+    def test_snapshot_density_matrix_pre_measure_det(self):
+        """Test snapshot density matrix before deterministic final measurement"""
+        shots = 10
+        label = "snap"
+        counts_targets = snapshot_state_counts_deterministic(shots)
+        statevec_targets = snapshot_state_pre_measure_statevector_deterministic(
+        )
+        circuits = snapshot_state_circuits_deterministic(label,
+                                                         'density_matrix',
+                                                         post_measure=False)
+
+        qobj = assemble(circuits, self.SIMULATOR, shots=shots)
+        job = self.SIMULATOR.run(qobj, backend_options=self.BACKEND_OPTS)
+        method = self.BACKEND_OPTS.get('method', 'automatic')
+        if method not in QasmSnapshotDensityMatrixTests.SUPPORTED_QASM_METHODS:
+            self.assertRaises(AerError, job.result)
+        else:
+            result = job.result()
+            self.is_completed(result)
+            self.compare_counts(result, circuits, counts_targets, delta=0)
+            # Check snapshots
+            for j, circuit in enumerate(circuits):
+                data = result.data(circuit)
+                snaps = self.density_snapshots(data, label)
+                self.assertTrue(len(snaps), 1)
+                target = np.outer(statevec_targets[j], statevec_targets[j].conj())
+                # Pre-measurement all memory bits should be 0
+                value = snaps.get('0x0')
+                self.assertTrue(np.allclose(value, target))
+
+    def test_snapshot_density_matrix_pre_measure_nondet(self):
+        """Test snapshot density matrix before non-deterministic final measurement"""
+        shots = 100
+        label = "snap"
+        counts_targets = snapshot_state_counts_nondeterministic(shots)
+        statevec_targets = snapshot_state_pre_measure_statevector_nondeterministic(
+        )
+        circuits = snapshot_state_circuits_nondeterministic(label,
+                                                            'density_matrix',
+                                                            post_measure=False)
+
+        qobj = assemble(circuits, self.SIMULATOR, shots=shots)
+        job = self.SIMULATOR.run(qobj, backend_options=self.BACKEND_OPTS)
+        method = self.BACKEND_OPTS.get('method', 'automatic')
+        if method not in QasmSnapshotDensityMatrixTests.SUPPORTED_QASM_METHODS:
+            self.assertRaises(AerError, job.result)
+        else:
+            result = job.result()
+            self.is_completed(result)
+            self.compare_counts(result,
+                                circuits,
+                                counts_targets,
+                                delta=0.2 * shots)
+            # Check snapshots
+            for j, circuit in enumerate(circuits):
+                data = result.data(circuit)
+                snaps = self.density_snapshots(data, label)
+                self.assertTrue(len(snaps), 1)
+                target = np.outer(statevec_targets[j], statevec_targets[j].conj())
+                value = snaps.get('0x0')
+                self.assertTrue(np.allclose(value, target))
+
+    def test_snapshot_density_matrix_post_measure_det(self):
+        """Test snapshot density matrix after deterministic final measurement"""
+        shots = 10
+        label = "snap"
+        counts_targets = snapshot_state_counts_deterministic(shots)
+        statevec_targets = snapshot_state_post_measure_statevector_deterministic(
+        )
+        circuits = snapshot_state_circuits_deterministic(label,
+                                                         'density_matrix',
+                                                         post_measure=True)
+
+        qobj = assemble(circuits, self.SIMULATOR, memory=True, shots=shots)
+        job = self.SIMULATOR.run(qobj, backend_options=self.BACKEND_OPTS)
+        method = self.BACKEND_OPTS.get('method', 'automatic')
+        if method not in QasmSnapshotDensityMatrixTests.SUPPORTED_QASM_METHODS:
+            self.assertRaises(AerError, job.result)
+        else:
+            result = job.result()
+            self.is_completed(result)
+            self.compare_counts(result, circuits, counts_targets, delta=0)
+            # Check snapshots
+            for i, circuit in enumerate(circuits):
+                data = result.data(circuit)
+                snaps = self.density_snapshots(data, label)
+                for j, mem in enumerate(data['memory']):
+                    target = statevec_targets[i].get(mem)
+                    target = np.outer(target, target.conj())
+                    value = snaps.get(mem)
+                    self.assertTrue(np.allclose(value, target))
+
+    def test_snapshot_density_matrix_post_measure_nondet(self):
+        """Test snapshot density matrix after non-deterministic final measurement"""
+        shots = 100
+        label = "snap"
+        counts_targets = snapshot_state_counts_nondeterministic(shots)
+        statevec_targets = snapshot_state_post_measure_statevector_nondeterministic(
+        )
+        circuits = snapshot_state_circuits_nondeterministic(label,
+                                                            'density_matrix',
+                                                            post_measure=True)
+
+        qobj = assemble(circuits, self.SIMULATOR, memory=True, shots=shots)
+        job = self.SIMULATOR.run(qobj, backend_options=self.BACKEND_OPTS)
+        method = self.BACKEND_OPTS.get('method', 'automatic')
+        if method not in QasmSnapshotDensityMatrixTests.SUPPORTED_QASM_METHODS:
+            self.assertRaises(AerError, job.result)
+        else:
+            result = job.result()
+            self.is_completed(result)
+            self.compare_counts(result,
+                                circuits,
+                                counts_targets,
+                                delta=0.2 * shots)
+            # Check snapshots
+            for i, circuit in enumerate(circuits):
+                data = result.data(circuit)
+                snaps = self.density_snapshots(data, label)
+                for j, mem in enumerate(data['memory']):
+                    target = statevec_targets[i].get(mem)
+                    target = np.outer(target, target.conj())
+                    value = snaps.get(mem)
+                    self.assertTrue(np.allclose(value, target))

test/terra/backends/test_qasm_density_matrix_simulator.py

@@ -46,6 +46,8 @@
 from test.terra.backends.qasm_simulator.qasm_method import QasmMethodTests
 # Snapshot tests
 from test.terra.backends.qasm_simulator.qasm_snapshot import QasmSnapshotStatevectorTests
+from test.terra.backends.qasm_simulator.qasm_snapshot import QasmSnapshotDensityMatrixTests
+from test.terra.backends.qasm_simulator.qasm_snapshot import QasmSnapshotStabilizerTests
 
 
 class TestQasmDensityMatrixSimulator(common.QiskitAerTestCase,
@@ -71,7 +73,9 @@ class TestQasmDensityMatrixSimulator(common.QiskitAerTestCase,
                                      QasmPauliNoiseTests,
                                      QasmResetNoiseTests,
                                      QasmKrausNoiseTests,
-                                     QasmSnapshotStatevectorTests):
+                                     QasmSnapshotStatevectorTests,
+                                     QasmSnapshotDensityMatrixTests,
+                                     QasmSnapshotStabilizerTests):
     """QasmSimulator density_matrix method tests."""
 
     BACKEND_OPTS = {