Feature/support density matrix simulation

docs/changelog.rst

@@ -1,6 +1,11 @@
 Changelog
 ~~~~~~~~~
 
+Unreleased
+--------------------
+
+* Add support for density matrix simulation.
+
 0.31.0 (October 2023)
 ---------------------
 

pytket/extensions/qulacs/backends/qulacs_backend.py

@@ -15,13 +15,13 @@
 """Methods to allow tket circuits to be ran on the Qulacs simulator
 """
 
-from typing import List, Optional, Sequence, Union, cast
+from typing import List, Optional, Sequence, Union, Type, cast
 from logging import warning
 from random import Random
 from uuid import uuid4
 import numpy as np
 from sympy import Expr
-from qulacs import Observable, QuantumState
+from qulacs import Observable, QuantumState, DensityMatrix
 from pytket.backends import (
     Backend,
     CircuitNotRunError,
@@ -106,7 +106,17 @@ class QulacsBackend(Backend):
         OpType.Barrier,
     }
 
-    def __init__(self) -> None:
+    def __init__(
+        self,
+        result_type: str = "state_vector",
+    ) -> None:
+        """
+        Backend for running simulations on the Qulacs simulator
+
+        :param result_type: Indicating the type of the simulation result
+            to be returned. It can be either "state_vector" or "density_matrix".
+            Defaults to "state_vector"
+        """
         super().__init__()
         self._backend_info = BackendInfo(
             type(self).__name__,
@@ -115,8 +125,16 @@ def __init__(self) -> None:
             Architecture([]),
             self._GATE_SET,
         )
-
-        self._sim = QuantumState
+        self._result_type = result_type
+        self._sim: Type[Union[QuantumState, DensityMatrix, "QuantumStateGpu"]]
+        if result_type == "state_vector":
+            self._sim = QuantumState
+        elif result_type == "density_matrix":
+            self._sim = DensityMatrix
+            self._supports_state = False
+            self._supports_density_matrix = True
+        else:
+            raise ValueError(f"Unsupported result type {result_type}")
 
     @property
     def _result_id_type(self) -> _ResultIdTuple:
@@ -193,7 +211,10 @@ def process_circuits(
                 circuit, reverse_index=True, replace_implicit_swaps=True
             )
             qulacs_circ.update_quantum_state(qulacs_state)
-            state = qulacs_state.get_vector()
+            if self._result_type == "state_vector":
+                state = qulacs_state.get_vector()  # type: ignore
+            else:
+                state = qulacs_state.get_matrix()  # type: ignore
             qubits = sorted(circuit.qubits, reverse=False)
             shots = None
             bits = None
@@ -214,28 +235,39 @@ def process_circuits(
                     )
                     shots = OutcomeArray.from_ints(samples, circuit.n_qubits)
                     shots = shots.choose_indices(choose_indices)
-            try:
-                phase = float(circuit.phase)
-                coeff = np.exp(phase * np.pi * 1j)
-                state *= coeff
-            except TypeError:
-                warning(
-                    "Global phase is dependent on a symbolic parameter, so cannot "
-                    "adjust for phase"
-                )
+            if self._result_type == "state_vector":
+                try:
+                    phase = float(circuit.phase)
+                    coeff = np.exp(phase * np.pi * 1j)
+                    state *= coeff
+                except TypeError:
+                    warning(
+                        "Global phase is dependent on a symbolic parameter, so cannot "
+                        "adjust for phase"
+                    )
             handle = ResultHandle(str(uuid4()))
-            self._cache[handle] = {
-                "result": BackendResult(
-                    state=state, shots=shots, c_bits=bits, q_bits=qubits
-                )
-            }
+            if self._result_type == "state_vector":
+                self._cache[handle] = {
+                    "result": BackendResult(
+                        state=state, shots=shots, c_bits=bits, q_bits=qubits
+                    )
+                }
+            else:
+                self._cache[handle] = {
+                    "result": BackendResult(
+                        density_matrix=state, shots=shots, c_bits=bits, q_bits=qubits
+                    )
+                }
             handle_list.append(handle)
             del qulacs_state
             del qulacs_circ
         return handle_list
 
     def _sample_quantum_state(
-        self, quantum_state: QuantumState, n_shots: int, rng: Optional[Random]
+        self,
+        quantum_state: Union[QuantumState, DensityMatrix, "QuantumStateGpu"],
+        n_shots: int,
+        rng: Optional[Random],
     ) -> List[int]:
         if rng:
             return quantum_state.sampling(n_shots, rng.randint(0, 2**32 - 1))
@@ -301,6 +333,9 @@ class QulacsGPUBackend(QulacsBackend):
         """
 
         def __init__(self) -> None:
+            """
+            Backend for running simulations on the Qulacs GPU simulator
+            """
             super().__init__()
             self._backend_info.name = type(self).__name__
             self._sim = QuantumStateGpu

tests/test_qulacs_backend.py

@@ -13,10 +13,11 @@
 # limitations under the License.
 
 from collections import Counter
-from typing import List, Sequence, Union, Optional
+from typing import List, Sequence, Union, Optional, Dict, Any
 import warnings
 import math
-from hypothesis import given, strategies
+from datetime import timedelta
+from hypothesis import given, strategies, settings
 import numpy as np
 import pytest
 from openfermion.ops import QubitOperator  # type: ignore
@@ -32,8 +33,10 @@
 
 def make_seeded_QulacsBackend(base: type[QulacsBackend]) -> type:
     class SeededQulacsBackend(base):  # type: ignore
-        def __init__(self, seed: int):
-            base.__init__(self)
+        def __init__(self, seed: int, kwargs: Optional[Dict[str, Any]] = None):
+            if kwargs is None:
+                kwargs = {}
+            base.__init__(self, **kwargs)
             self._seed = seed
 
         def process_circuits(
@@ -50,7 +53,12 @@ def process_circuits(
     return SeededQulacsBackend
 
 
-backends = [QulacsBackend(), make_seeded_QulacsBackend(QulacsBackend)(-1)]
+backends = [
+    QulacsBackend(),
+    make_seeded_QulacsBackend(QulacsBackend)(-1),
+    QulacsBackend(result_type="density_matrix"),
+    make_seeded_QulacsBackend(QulacsBackend)(-1, {"result_type": "density_matrix"}),
+]
 
 try:
     from pytket.extensions.qulacs import QulacsGPUBackend
@@ -100,6 +108,15 @@ def h2_4q_circ(theta: float) -> Circuit:
     return circ
 
 
+def test_properties() -> None:
+    svb = QulacsBackend()
+    dmb = QulacsBackend(result_type="density_matrix")
+    assert not svb.supports_density_matrix
+    assert svb.supports_state
+    assert not dmb.supports_state
+    assert dmb.supports_density_matrix
+
+
 def test_get_state() -> None:
     qulacs_circ = h2_4q_circ(PARAM)
     correct_state = np.array(
@@ -124,12 +141,20 @@ def test_get_state() -> None:
     )
     for b in backends:
         qulacs_circ = b.get_compiled_circuit(qulacs_circ)
-        qulacs_state = b.run_circuit(qulacs_circ).get_state()
-        assert np.allclose(qulacs_state, correct_state)
+        if b.supports_state:
+            qulacs_state = b.run_circuit(qulacs_circ).get_state()
+            assert np.allclose(qulacs_state, correct_state)
+        if b.supports_density_matrix:
+            qulacs_state = b.run_circuit(qulacs_circ).get_density_matrix()
+            assert np.allclose(
+                qulacs_state, np.outer(correct_state, correct_state.conj())
+            )
 
 
 def test_statevector_phase() -> None:
     for b in backends:
+        if not b.supports_state:
+            continue
         circ = Circuit(2)
         circ.H(0).CX(0, 1)
         circ = b.get_compiled_circuit(circ)
@@ -151,14 +176,24 @@ def test_swaps_basisorder() -> None:
         assert c.n_gates_of_type(OpType.CX) == 1
         c = b.get_compiled_circuit(c)
         res = b.run_circuit(c)
-        s_ilo = res.get_state(basis=BasisOrder.ilo)
-        s_dlo = res.get_state(basis=BasisOrder.dlo)
-        correct_ilo = np.zeros((16,))
-        correct_ilo[4] = 1.0
-        assert np.allclose(s_ilo, correct_ilo)
-        correct_dlo = np.zeros((16,))
-        correct_dlo[2] = 1.0
-        assert np.allclose(s_dlo, correct_dlo)
+        if b.supports_state:
+            s_ilo = res.get_state(basis=BasisOrder.ilo)
+            s_dlo = res.get_state(basis=BasisOrder.dlo)
+            correct_ilo = np.zeros((16,))
+            correct_ilo[4] = 1.0
+            assert np.allclose(s_ilo, correct_ilo)
+            correct_dlo = np.zeros((16,))
+            correct_dlo[2] = 1.0
+            assert np.allclose(s_dlo, correct_dlo)
+        if b.supports_density_matrix:
+            s_ilo = res.get_density_matrix(basis=BasisOrder.ilo)
+            s_dlo = res.get_density_matrix(basis=BasisOrder.dlo)
+            correct_ilo = np.zeros((16,))
+            correct_ilo[4] = 1.0
+            assert np.allclose(s_ilo, np.outer(correct_ilo, correct_ilo.conj()))
+            correct_dlo = np.zeros((16,))
+            correct_dlo[2] = 1.0
+            assert np.allclose(s_dlo, np.outer(correct_dlo, correct_dlo.conj()))
 
 
 def from_OpenFermion(openf_op: "QubitOperator") -> QubitPauliOperator:
@@ -204,8 +239,18 @@ def test_basisorder() -> None:
         c.X(1)
         b.process_circuit(c)
         res = b.run_circuit(c)
-        assert (res.get_state() == np.asarray([0, 1, 0, 0])).all()
-        assert (res.get_state(basis=BasisOrder.dlo) == np.asarray([0, 0, 1, 0])).all()
+        if b.supports_state:
+            assert (res.get_state() == np.asarray([0, 1, 0, 0])).all()
+            assert (
+                res.get_state(basis=BasisOrder.dlo) == np.asarray([0, 0, 1, 0])
+            ).all()
+        if b.supports_density_matrix:
+            sv = np.asarray([0, 1, 0, 0])
+            assert (res.get_density_matrix() == np.outer(sv, sv.conj())).all()
+            sv = np.asarray([0, 0, 1, 0])
+            assert (
+                res.get_density_matrix(basis=BasisOrder.dlo) == np.outer(sv, sv.conj())
+            ).all()
         c.measure_all()
         res = b.run_circuit(c, n_shots=4, seed=4)
         assert res.get_shots().shape == (4, 2)
@@ -272,22 +317,29 @@ def test_backend_with_circuit_permutation() -> None:
     for b in backends:
         c = Circuit(3).X(0).SWAP(0, 1).SWAP(0, 2)
         qubits = c.qubits
-        sv = b.run_circuit(c).get_state()
+        if b.supports_state:
+            sv = b.run_circuit(c).get_state()
+        else:
+            sv = b.run_circuit(c).get_density_matrix()
         # convert swaps to implicit permutation
         c.replace_SWAPs()
         assert c.implicit_qubit_permutation() == {
             qubits[0]: qubits[1],
             qubits[1]: qubits[2],
             qubits[2]: qubits[0],
         }
-        sv1 = b.run_circuit(c).get_state()
+        if b.supports_state:
+            sv1 = b.run_circuit(c).get_state()
+        else:
+            sv1 = b.run_circuit(c).get_density_matrix()
         assert np.allclose(sv, sv1, atol=1e-10)
 
 
 @given(
     n_shots=strategies.integers(min_value=1, max_value=10),  # type: ignore
     n_bits=strategies.integers(min_value=0, max_value=10),
 )
+@settings(deadline=timedelta(seconds=1))
 def test_shots_bits_edgecases(n_shots, n_bits) -> None:
     c = Circuit(n_bits, n_bits)