diff --git a/docs/changelog.rst b/docs/changelog.rst index 2b42b12..6b9fd9b 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,6 +1,11 @@ Changelog ~~~~~~~~~ +Unreleased +-------------------- + +* Add support for density matrix simulation. + 0.31.0 (October 2023) --------------------- diff --git a/pytket/extensions/qulacs/backends/qulacs_backend.py b/pytket/extensions/qulacs/backends/qulacs_backend.py index 1a7385e..dfbc77c 100644 --- a/pytket/extensions/qulacs/backends/qulacs_backend.py +++ b/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 diff --git a/tests/test_qulacs_backend.py b/tests/test_qulacs_backend.py index 9ccaf7f..bcf72c7 100644 --- a/tests/test_qulacs_backend.py +++ b/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,7 +317,10 @@ 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() == { @@ -280,7 +328,10 @@ def test_backend_with_circuit_permutation() -> None: 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) @@ -288,6 +339,7 @@ def test_backend_with_circuit_permutation() -> None: 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)