Include density_matrix flag in random_clifford and random_pauli

src/qibo/quantum_info/random_ensembles.py

@@ -557,6 +557,7 @@ def random_density_matrix(
 def random_clifford(
     nqubits: int,
     return_circuit: bool = True,
+    density_matrix: bool = False,
     seed=None,
     backend=None,
 ):
@@ -567,6 +568,8 @@ def random_clifford(
         nqubits (int): number of qubits.
         return_circuit (bool, optional): if ``True``, returns a :class:`qibo.models.Circuit`
             object. If ``False``, returns an ``ndarray`` object. Defaults to ``False``.
+        density_matrix (bool, optional): used when ``return_circuit=True``. If `True`,
+            the circuit would evolve density matrices. Defaults to ``False``.
         seed (int or :class:`numpy.random.Generator`, optional): Either a generator of
             random numbers or a fixed seed to initialize a generator. If ``None``,
             initializes a generator with a random seed. Defaults to ``None``.
@@ -672,7 +675,9 @@ def random_clifford(
                 delta_matrix[k, j] = b
 
     # get first element of the Borel group
-    clifford_circuit = _operator_from_hadamard_free_group(gamma_matrix, delta_matrix)
+    clifford_circuit = _operator_from_hadamard_free_group(
+        gamma_matrix, delta_matrix, density_matrix
+    )
 
     # Apply permutated Hadamard layer
     for qubit, had in enumerate(hadamards):
@@ -683,7 +688,15 @@ def random_clifford(
     clifford_circuit += _operator_from_hadamard_free_group(
         gamma_matrix_prime,
         delta_matrix_prime,
-        random_pauli(nqubits, depth=1, return_circuit=True, seed=seed, backend=backend),
+        density_matrix,
+        random_pauli(
+            nqubits,
+            depth=1,
+            return_circuit=True,
+            density_matrix=density_matrix,
+            seed=seed,
+            backend=backend,
+        ),
     )
 
     if return_circuit is False:
@@ -698,6 +711,7 @@ def random_pauli(
     max_qubits: Optional[int] = None,
     subset: Optional[list] = None,
     return_circuit: bool = True,
+    density_matrix: bool = False,
     seed=None,
     backend=None,
 ):
@@ -720,6 +734,8 @@ def random_pauli(
         return_circuit (bool, optional): if ``True``, returns a :class:`qibo.models.Circuit`
             object. If ``False``, returns an ``ndarray`` with shape (qubits, depth, 2, 2)
             that contains all Pauli matrices that were sampled. Defaults to ``True``.
+        density_matrix (bool, optional): used when ``return_circuit=True``. If `True`,
+            the circuit would evolve density matrices. Defaults to ``False``.
         seed (int or :class:`numpy.random.Generator`, optional): Either a generator of
             random numbers or a fixed seed to initialize a generator. If ``None``,
             initializes a generator with a random seed. Defaults to ``None``.
@@ -819,7 +835,7 @@ def random_pauli(
     indexes = [[keys[item] for item in row] for row in indexes]
 
     if return_circuit:
-        gate_grid = Circuit(max_qubits)
+        gate_grid = Circuit(max_qubits, density_matrix=density_matrix)
         for qubit, row in zip(qubits, indexes):
             for column_item in row:
                 if subset[column_item] != gates.I:
@@ -1100,14 +1116,18 @@ def _sample_from_quantum_mallows_distribution(nqubits: int, local_state):
     return hadamards, permutations
 
 
-def _operator_from_hadamard_free_group(gamma_matrix, delta_matrix, pauli_operator=None):
+def _operator_from_hadamard_free_group(
+    gamma_matrix, delta_matrix, density_matrix: bool = False, pauli_operator=None
+):
     """Calculates an element :math:`F` of the Hadamard-free group :math:`\\mathcal{F}_{n}`,
     where :math:`n` is the number of qubits ``nqubits``. For more details,
     see Reference [1].
 
     Args:
         gamma_matrix (ndarray): :math:`\\, n \\times n \\,` binary matrix.
         delta_matrix (ndarray): :math:`\\, n \\times n \\,` binary matrix.
+        density_matrix (bool, optional): used when ``return_circuit=True``. If `True`,
+            the circuit would evolve density matrices. Defaults to ``False``.
         pauli_operator (:class:`qibo.models.Circuit`, optional): a :math:`n`-qubit
             Pauli operator. If ``None``, it is assumed to be the Identity.
             Defaults to ``None``.
@@ -1128,7 +1148,7 @@ def _operator_from_hadamard_free_group(gamma_matrix, delta_matrix, pauli_operato
         )
 
     nqubits = len(gamma_matrix)
-    circuit = Circuit(nqubits)
+    circuit = Circuit(nqubits, density_matrix=density_matrix)
 
     if pauli_operator is not None:
         circuit += pauli_operator

tests/test_quantum_info_random.py

@@ -273,9 +273,10 @@ def test_random_density_matrix(backend, dims, pure, metric, basis, normalize):
 
 
 @pytest.mark.parametrize("seed", [10])
+@pytest.mark.parametrize("density_matrix", [False, True])
 @pytest.mark.parametrize("return_circuit", [True, False])
 @pytest.mark.parametrize("nqubits", [1, 2])
-def test_random_clifford(backend, nqubits, return_circuit, seed):
+def test_random_clifford(backend, nqubits, return_circuit, density_matrix, seed):
     with pytest.raises(TypeError):
         test = random_clifford(
             nqubits="1", return_circuit=return_circuit, backend=backend
@@ -301,7 +302,11 @@ def test_random_clifford(backend, nqubits, return_circuit, seed):
     result = backend.cast(result, dtype=result.dtype)
 
     matrix = random_clifford(
-        nqubits, return_circuit=return_circuit, seed=seed, backend=backend
+        nqubits,
+        return_circuit=return_circuit,
+        density_matrix=density_matrix,
+        seed=seed,
+        backend=backend,
     )
 
     if return_circuit:
@@ -374,8 +379,11 @@ def test_pauli_single(backend):
 @pytest.mark.parametrize("max_qubits", [None])
 @pytest.mark.parametrize("subset", [None, ["I", "X"]])
 @pytest.mark.parametrize("return_circuit", [True, False])
+@pytest.mark.parametrize("density_matrix", [False, True])
 @pytest.mark.parametrize("seed", [10])
-def test_random_pauli(backend, qubits, depth, max_qubits, subset, return_circuit, seed):
+def test_random_pauli(
+    backend, qubits, depth, max_qubits, subset, return_circuit, density_matrix, seed
+):
     result_complete_set = np.array(
         [
             [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
@@ -390,7 +398,7 @@ def test_random_pauli(backend, qubits, depth, max_qubits, subset, return_circuit
     result_subset = backend.identity_density_matrix(2, normalize=False)
 
     matrix = random_pauli(
-        qubits, depth, max_qubits, subset, return_circuit, seed, backend
+        qubits, depth, max_qubits, subset, return_circuit, density_matrix, seed, backend
     )
 
     if return_circuit: