Use the atol option in qutip to determine the tolerance for measureme…

…nt. And filter out the states with a probability smaller than this.

src/qutip_qip/operations/measurement.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 
+import qutip
 from qutip import basis, Qobj
 from qutip.measurement import measurement_statistics
 from .gates import expand_operator
@@ -77,7 +78,6 @@ def measurement_comp_basis(self, state):
 
         n = int(np.log2(state.shape[0]))
         target = self.targets[0]
-
         if target < n:
             op0 = basis(2, 0) * basis(2, 0).dag()
             op1 = basis(2, 1) * basis(2, 1).dag()
@@ -89,7 +89,22 @@ def measurement_comp_basis(self, state):
             expand_operator(op, N=n, targets=self.targets)
             for op in measurement_ops
         ]
-        return measurement_statistics(state, measurement_ops)
+
+        try:
+            # qutip-v5
+            measurement_tol = qutip.settings.core["atol"] ** 2
+        except AttributeError:
+            # qutip-v4
+            measurement_tol = qutip.settings.atol**2
+        states, probabilities = measurement_statistics(state, measurement_ops)
+        probabilities = [
+            p if p > measurement_tol else 0.0 for p in probabilities
+        ]
+        states = [
+            s if p > measurement_tol else None
+            for s, p in zip(states, probabilities)
+        ]
+        return states, probabilities
 
     def __str__(self):
         str_name = ("Measurement(%s, target=%s, classical_store=%s)") % (

tests/test_measurement.py

@@ -8,6 +8,7 @@
 from qutip.measurement import (measure_povm, measurement_statistics_povm,
                                 measure_observable,
                                 measurement_statistics_observable)
+import qutip
 
 
 @pytest.mark.repeat(10)
@@ -67,3 +68,11 @@ def test_measurement_collapse(index):
             states_11, probability_11 = Mprime.measurement_comp_basis(state)
             assert probability_11[1] == 1
             assert states_11[0] is None
+
+
+def test_against_numerical_error():
+    state = qutip.Qobj([[1], [1.e-12]])
+    measurement = Measurement("M", 0)
+    states, probabilites = measurement.measurement_comp_basis(state)
+    assert states[1] is None
+    assert probabilites[1] == 0.