Repository: qiboteam/qibo

Issue:
ValueError when probability vector sums to 1 up to machine precision
```
from itertools import product
from qibo import gates
from qibo.models import Circuit
from qibo.hamiltonians import SymbolicHamiltonian
from qibo.symbols import I, X, Y, Z


def get_pauli_basis(qubits: int) -> list:
    ''' Returns list with all 4^qubits Pauli observables '''
    pauli_single = [I, X, Y, Z]
    aux = list(product(range(len(pauli_single)), repeat = qubits))
    observables = list()
    for row in aux:
        observable = 1
        for qubit, pauli in enumerate(row):
            observable = observable * pauli_single[pauli](qubit)
        observables.append(observable)
    return [SymbolicHamiltonian(pauli) for pauli in observables]


def get_probabilities(qubits: int, param: float = 1.) -> list:
    ''' Returns probability vector for a generalized depolarizing channel '''
    d = 2**qubits
    num_terms = d**2
    max_param = num_terms / (num_terms - 1)
    prob_identity = 1 - param / max_param
    prob_pauli = param / num_terms
    return [prob_identity] + [prob_pauli] * (num_terms - 1)


qubits = 2
param = 0.006

paulis = get_pauli_basis(qubits)
paulis = [pauli.dense.matrix for pauli in paulis]
paulis = list(product([(0, 1)], paulis))

probs = get_probabilities(qubits, param)

circuit = Circuit(qubits)
circuit.add((gates.UnitaryChannel(probs, paulis)))
```


This minimal code example raises the following:

```ValueError: UnitaryChannel probability sum should be between 0 and 1 but is 1.0000000000000002.```

Discussion:

Patch:
diff --git a/src/qibo/abstractions/gates.py b/src/qibo/abstractions/gates.py
index caf8033d5..cd7d35120 100644
--- a/src/qibo/abstractions/gates.py
+++ b/src/qibo/abstractions/gates.py
@@ -1467,10 +1467,6 @@ class UnitaryChannel(KrausChannel):
         self.name = "UnitaryChannel"
         self.probs = p
         self.psum = sum(p)
-        if self.psum > 1 or self.psum <= 0:
-            raise_error(ValueError, "UnitaryChannel probability sum should be "
-                                    "between 0 and 1 but is {}."
-                                    "".format(self.psum))
         self.seed = seed
         self.density_matrix = False
         self.init_args = [p, self.gates]
diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py
index 5c5c3a7d4..5c51177c9 100644
--- a/src/qibo/backends/abstract.py
+++ b/src/qibo/backends/abstract.py
@@ -15,6 +15,8 @@ class AbstractBackend(ABC):
         self.available_platforms = [None]
 
         self.precision = 'double'
+        from qibo.config import PRECISION_TOL_DOUBLE
+        self.precision_tol = PRECISION_TOL_DOUBLE
         self._dtypes = {'DTYPEINT': 'int64', 'DTYPE': 'float64',
                         'DTYPECPX': 'complex128'}
 
@@ -56,11 +58,15 @@ class AbstractBackend(ABC):
 
     def set_precision(self, dtype):
         if dtype == 'single':
+            from qibo.config import PRECISION_TOL_SINGLE
             self._dtypes['DTYPE'] = 'float32'
             self._dtypes['DTYPECPX'] = 'complex64'
+            self.precision_tol = PRECISION_TOL_SINGLE
         elif dtype == 'double':
+            from qibo.config import PRECISION_TOL_DOUBLE
             self._dtypes['DTYPE'] = 'float64'
             self._dtypes['DTYPECPX'] = 'complex128'
+            self.precision_tol = PRECISION_TOL_DOUBLE
         else:
             raise_error(ValueError, f'dtype {dtype} not supported.')
         self.precision = dtype
diff --git a/src/qibo/config.py b/src/qibo/config.py
index 81cbdc907..6644c55db 100644
--- a/src/qibo/config.py
+++ b/src/qibo/config.py
@@ -21,6 +21,10 @@ EINSUM_CHARS = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
 # Eigenvalues smaller than this cut-off are ignored in entropy calculation
 EIGVAL_CUTOFF = 1e-14
 
+# Tolerance for the probability sum check in the unitary channel
+PRECISION_TOL_SINGLE = 1e-8
+PRECISION_TOL_DOUBLE = 1e-12
+
 # Batch size for sampling shots in measurement frequencies calculation
 SHOT_BATCH_SIZE = 2 ** 18
 
diff --git a/src/qibo/core/gates.py b/src/qibo/core/gates.py
index 05ed72557..b5b5bb14f 100644
--- a/src/qibo/core/gates.py
+++ b/src/qibo/core/gates.py
@@ -927,6 +927,10 @@ class UnitaryChannel(KrausChannel, abstract_gates.UnitaryChannel):
                  seed: Optional[int] = None):
         BackendGate.__init__(self)
         abstract_gates.UnitaryChannel.__init__(self, p, ops, seed=seed)
+        if self.psum > 1 + K.precision_tol or self.psum <= 0:
+            raise_error(ValueError, "UnitaryChannel probability sum should be "
+                                    "between 0 and 1 but is {}."
+                                    "".format(self.psum))
         self.set_seed()
 
     def calculate_inverse_gates(self):


Test patch:
diff --git a/src/qibo/tests/test_abstract_gates.py b/src/qibo/tests/test_abstract_gates.py
index 363aa1a3f..525ab95c6 100644
--- a/src/qibo/tests/test_abstract_gates.py
+++ b/src/qibo/tests/test_abstract_gates.py
@@ -247,8 +247,6 @@ def test_unitary_channel_init():
         gate = gates.UnitaryChannel(2 * [0.1], ops)
     with pytest.raises(ValueError):
         gate = gates.UnitaryChannel(4 * [-0.1], ops)
-    with pytest.raises(ValueError):
-        gate = gates.UnitaryChannel(4 * [0.5], ops)
 
 
 def test_pauli_noise_channel_init():
diff --git a/src/qibo/tests/test_backends_init.py b/src/qibo/tests/test_backends_init.py
index 7a00673f4..6837fd199 100644
--- a/src/qibo/tests/test_backends_init.py
+++ b/src/qibo/tests/test_backends_init.py
@@ -52,10 +52,13 @@ def test_set_precision(backend, precision):
     backends.set_precision(precision)
     if precision == "single":
         expected_dtype = K.backend.complex64
+        expected_tol = 1e-8
     else:
         expected_dtype = K.backend.complex128
+        expected_tol = 1e-12
     assert backends.get_precision() == precision
     assert K.dtypes('DTYPECPX') == expected_dtype
+    assert K.precision_tol == expected_tol
     # Test that circuits use proper precision
     circuit = models.Circuit(2)
     circuit.add([gates.H(0), gates.H(1)])
diff --git a/src/qibo/tests/test_core_gates.py b/src/qibo/tests/test_core_gates.py
index 6b41c17e4..14174f9cc 100644
--- a/src/qibo/tests/test_core_gates.py
+++ b/src/qibo/tests/test_core_gates.py
@@ -423,7 +423,7 @@ def test_general_channel(backend):
     initial_rho = random_density_matrix(2)
     gate = gates.KrausChannel([((1,), a1), ((0, 1), a2)])
     assert gate.target_qubits == (0, 1)
-    final_rho = gate(np.copy(initial_rho))
+    final_rho = gate(K.cast(np.copy(initial_rho)))
     m1 = np.kron(np.eye(2), K.to_numpy(a1))
     m2 = K.to_numpy(a2)
     target_rho = (m1.dot(initial_rho).dot(m1.conj().T) +
@@ -476,6 +476,28 @@ def test_unitary_channel(backend):
     K.assert_allclose(final_state, target_state)
 
 
+@pytest.mark.parametrize("precision", ["double", "single"])
+def test_unitary_channel_probability_tolerance(backend, precision):
+    """Create ``UnitaryChannel`` with probability sum within tolerance (see #562)."""
+    import qibo
+    original_precision = qibo.get_precision()
+    qibo.set_precision(precision)
+    nqubits = 2
+    param = 0.006
+    num_terms = 2 ** (2 * nqubits)
+    max_param = num_terms / (num_terms - 1)
+    prob_identity = 1 - param / max_param
+    prob_pauli = param / num_terms
+    probs = [prob_identity] + [prob_pauli] * (num_terms - 1)
+    if precision == "double":
+        probs = np.array(probs, dtype="float64")
+    else:
+        probs = np.array(probs, dtype="float32")
+    matrices = len(probs) * [((0, 1), np.random.random((4, 4)))]
+    gate = gates.UnitaryChannel(probs, matrices)
+    qibo.set_precision(original_precision)
+
+
 def test_unitary_channel_errors():
     """Check errors raised by ``gates.UnitaryChannel``."""
     a1 = np.array([[0, 1], [1, 0]])
@@ -488,7 +510,7 @@ def test_unitary_channel_errors():
     # Probability > 1
     with pytest.raises(ValueError):
         gate = gates.UnitaryChannel([1.1, 0.2], matrices)
-    # Probability sum = 0
+    # Probability sum < 0
     with pytest.raises(ValueError):
         gate = gates.UnitaryChannel([0.0, 0.0], matrices)