Merge pull request #1272 from qiboteam/qulacs

`QulacsBackend`

Author: BrunoLiegiBastonLiegi

doc/source/api-reference/qibo.rst

@@ -2440,6 +2440,16 @@ The default backend order is qibojit (if available), tensorflow (if available),
 numpy. The default backend can be changed using the ``QIBO_BACKEND`` environment
 variable.
 
+Qibo optionally provides an interface to `qulacs <https://github.com/qulacs/qulacs>`_ through the :class:`qibo.backends.qulacs.QulacsBackend`. To use ``qulacs`` for simulating a quantum circuit you can globally set the backend as in the other cases
+
+.. testcode:: python
+
+   import qibo
+   qibo.set_backend("qulacs")
+
+.. note::
+   GPU simulation through ``qulacs`` is not supported yet.
+
 .. autoclass:: qibo.backends.abstract.Backend
     :members:
     :member-order: bysource

doc/source/getting-started/backends.png

@@ -40,6 +40,7 @@ if the corresponding packages are installed, following the hierarchy below:
 * :ref:`installing-pytorch`: a pure PyTorch implementation for quantum simulation which provides access to gradient descent optimization and the possibility to implement classical and quantum architectures together. This backend is not optimized for memory and speed, use :ref:`installing-qibojit` instead.
 * :ref:`clifford <Clifford>`: a specialized backend for the simulation of quantum circuits with Clifford gates. This backend uses :ref:`installing-qibojit` and/or :ref:`installing-numpy`.
 * `qibotn <https://qibo.science/qibotn/stable/>`_: an interface to Tensor Networks simulation algorithms designed for GPUs and multi-node CPUs. This backend makes possible scaling quantum circuit simulation to a larger number of qubits.
+* `qulacs <https://github.com/qulacs/qulacs>`_: an interface to the `qulacs` library for quantum simulation. GPU support is not available yet.
 
 The default backend that is used is the first available from the above list.
 The user can switch to a different using the ``qibo.set_backend`` method

doc/source/getting-started/backends.rst

@@ -35,6 +35,7 @@ numpy = "^1.26.4"
 networkx = "^3.2.1"
 tensorflow = { version = "^2.16.1", markers = "sys_platform == 'linux' or sys_platform == 'darwin'", optional = true }
 torch = { version = "^2.1.1", optional = true }
+qulacs = { version = "^0.6.4", optional = true, markers="(sys_platform == 'darwin' and python_version > '3.9') or sys_platform != 'darwin'"}
 
 [tool.poetry.group.dev]
 optional = true
@@ -55,6 +56,7 @@ sphinx-markdown-tables = "^0.0.17"
 sphinx-copybutton = "^0.5.2"
 nbsphinx = "^0.8.12"
 ipython = "^8.10.0"
+qulacs = { version = "^0.6.4", markers="(sys_platform == 'darwin' and python_version > '3.9') or sys_platform != 'darwin'"}
 seaborn = "^0.13.2"
 ipykernel = "^6.29.4"
 
@@ -74,6 +76,7 @@ torch = "^2.1.1"
 qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 stim = "^1.12.0"
+qulacs = { version = "^0.6.4", markers="(sys_platform == 'darwin' and python_version > '3.9') or sys_platform != 'darwin'" }
 
 [tool.poe.tasks]
 test = "pytest"

doc/source/qibo_ecosystem.png

@@ -11,7 +11,7 @@
 from qibo.backends.tensorflow import TensorflowBackend
 from qibo.config import log, raise_error
 
-QIBO_NATIVE_BACKENDS = ("numpy", "tensorflow", "pytorch")
+QIBO_NATIVE_BACKENDS = ("numpy", "tensorflow", "pytorch", "qulacs")
 QIBO_NON_NATIVE_BACKENDS = ("qibojit", "qibolab", "qibo-cloud-backends", "qibotn")
 
 
@@ -39,6 +39,10 @@ def load(backend: str, **kwargs) -> Backend:
             engine = kwargs.pop("platform", None)
             kwargs["engine"] = engine
             return CliffordBackend(**kwargs)
+        elif backend == "qulacs":
+            from qibo.backends.qulacs import QulacsBackend
+
+            return QulacsBackend()
         else:
             raise_error(
                 ValueError,

poetry.lock

@@ -0,0 +1,86 @@
+import re
+
+import numpy as np
+import qulacs  # pylint: disable=import-error
+from qulacs import (  # pylint: disable=no-name-in-module, import-error
+    QuantumCircuitSimulator,
+    converter,
+)
+
+from qibo import __version__
+from qibo.backends import NumpyBackend
+from qibo.config import raise_error
+from qibo.result import CircuitResult, QuantumState
+
+
+def circuit_to_qulacs(
+    circuit: "qibo.Circuit",
+) -> "qulacs.QuantumCircuit":  # pylint: disable=no-member
+    """
+    Converts a qibo circuit in a qulacs circuit.
+
+    Args:
+        circuit (:class:`qibo.models.circuit.Circuit`): Input circuit to convert.
+
+    Returns:
+        qulacs.QuantumCircuit: The converted qulacs circuit.
+    """
+    qasm_str = re.sub("^//.+\n", "", circuit.to_qasm())
+    qasm_str = re.sub(r"creg\s.+;", "", qasm_str)
+    qasm_str = re.sub(r"measure\s.+;", "", qasm_str)
+    circ = converter.convert_QASM_to_qulacs_circuit(qasm_str.splitlines())
+    return circ
+
+
+class QulacsBackend(NumpyBackend):
+
+    def __init__(self):
+        super().__init__()
+
+        self.name = "qulacs"
+        self.versions = {"qibo": __version__, "qulacs": qulacs.__version__}
+        self.device = "CPU"
+
+    def execute_circuit(
+        self,
+        circuit: "qibo.Circuit",
+        initial_state=None,
+        nshots: int = 1000,
+    ):
+        """Execute a circuit with qulacs.
+
+        Args:
+            circuit (:class:`qibo.models.circuit.Circuit`): Input circuit.
+            nshots (int, optional): Number of shots to perform if ``circuit`` has measurements.
+                Defaults to :math:`10^{3}`.
+
+        Returns:
+            :class:`qibo.result.CircuitResult`: Object storing to the final results.
+        """
+        if initial_state is not None:
+            raise_error(
+                NotImplementedError,
+                "The use of an initial state is not supported yet by the `QulacsBackend`.",
+            )
+        circ = circuit_to_qulacs(circuit)
+        state = (
+            qulacs.DensityMatrix(circuit.nqubits)  # pylint: disable=no-member
+            if circuit.density_matrix
+            else qulacs.QuantumState(circuit.nqubits)  # pylint: disable=no-member
+        )
+        sim = QuantumCircuitSimulator(circ, state)
+        sim.simulate()
+        if circuit.density_matrix:
+            dim = 2**circuit.nqubits
+            state = (
+                state.get_matrix()
+                .reshape(2 * circuit.nqubits * (2,))
+                .T.reshape(dim, dim)
+            )
+        else:
+            state = state.get_vector().reshape(circuit.nqubits * (2,)).T.ravel()
+        if len(circuit.measurements) > 0:
+            return CircuitResult(
+                state, circuit.measurements, backend=self, nshots=nshots
+            )
+        return QuantumState(state, backend=self)

pyproject.toml

@@ -322,7 +322,9 @@ def execute(self, backend=None):
 
             # Calling and executing the quantum circuit
             qc = self.construct_qae_circuit(k_i)
-            samples = qc(nshots=n_shots_i).frequencies(binary=True)["1"]
+            samples = backend.execute_circuit(qc, nshots=n_shots_i).frequencies(
+                binary=True
+            )["1"]
 
             samples_history.append(samples)
             n_shots_history.append(n_shots_i)

src/qibo/backends/__init__.py

@@ -46,6 +46,12 @@ def get_backend(backend_name):
     except (ModuleNotFoundError, ImportError):
         pass
 
+try:
+    get_backend("qulacs")
+    QULACS_INSTALLED = True
+except ModuleNotFoundError:
+    QULACS_INSTALLED = False
+
 
 def pytest_runtest_setup(item):
     ALL = {"darwin", "linux"}
@@ -54,6 +60,9 @@ def pytest_runtest_setup(item):
     if supported_platforms and plat not in supported_platforms:  # pragma: no cover
         # case not covered by workflows
         pytest.skip(f"Cannot run test on platform {plat}.")
+    elif not QULACS_INSTALLED and item.fspath.purebasename == "test_backends_qulacs":
+        # case not covered by workflows
+        pytest.skip(f"Cannot test `qulacs` on platform {plat}.")
 
 
 def pytest_configure(config):

src/qibo/backends/qulacs.py

@@ -1,4 +1,5 @@
 import platform
+import sys
 
 import numpy as np
 import pytest
@@ -113,10 +114,14 @@ def test_construct_backend(backend):
 
 def test_list_available_backends():
     tensorflow = False if platform.system() == "Windows" else True
+    qulacs = (
+        False if platform.system() == "Darwin" and sys.version_info[1] == 9 else True
+    )
     available_backends = {
         "numpy": True,
         "tensorflow": tensorflow,
         "pytorch": True,
+        "qulacs": qulacs,
         "qibojit": {"numba": True, "cupy": False, "cuquantum": False},
         "qibolab": False,
         "qibo-cloud-backends": False,

src/qibo/models/iqae.py

@@ -0,0 +1,43 @@
+import random
+
+import numpy as np
+import pytest
+
+from qibo import Circuit, gates
+from qibo.backends import GlobalBackend, MetaBackend, NumpyBackend, set_backend
+from qibo.quantum_info import random_clifford, random_density_matrix, random_statevector
+
+numpy_bkd = NumpyBackend()
+
+
+@pytest.mark.parametrize("density_matrix", [True, False])
+@pytest.mark.parametrize("with_measurements", [True, False])
+def test_qulacs(density_matrix, with_measurements):
+    c = random_clifford(3, backend=numpy_bkd, density_matrix=density_matrix)
+    if with_measurements:
+        measured_qubits = random.sample([0, 1, 2], 2)
+        c.add(gates.M(*measured_qubits))
+    qulacs_bkd = MetaBackend.load("qulacs")
+    nshots = 1000
+    qulacs_res = qulacs_bkd.execute_circuit(c, nshots=nshots)
+    numpy_res = numpy_bkd.execute_circuit(c, nshots=nshots)
+    numpy_bkd.assert_allclose(numpy_res.probabilities(), qulacs_res.probabilities())
+    if with_measurements:
+        numpy_freq = numpy_res.frequencies(binary=True)
+        qulacs_freq = qulacs_res.frequencies(binary=True)
+        numpy_freq = [numpy_freq.get(state, 0) / nshots for state in range(8)]
+        qulacs_freq = [qulacs_freq.get(state, 0) / nshots for state in range(8)]
+        numpy_bkd.assert_allclose(numpy_freq, qulacs_freq, atol=1e-1)
+
+
+def test_initial_state_error():
+    c = Circuit(1)
+    qulacs_bkd = MetaBackend.load("qulacs")
+    initial_state = np.array([0.0, 1.0])
+    with pytest.raises(NotImplementedError):
+        qulacs_bkd.execute_circuit(c, initial_state=initial_state)
+
+
+def test_set_backend():
+    set_backend("qulacs")
+    assert GlobalBackend().name == "qulacs"

tests/conftest.py

@@ -424,7 +424,7 @@ def test_measurement_basis(backend, nqubits, outcome):
         c.add(gates.X(q) for q in range(nqubits))
     c.add(gates.H(q) for q in range(nqubits))
     c.add(gates.M(*range(nqubits), basis=gates.X))
-    result = c(nshots=100)
+    result = backend.execute_circuit(c, nshots=100)
     assert result.frequencies() == {nqubits * str(outcome): 100}
 
 
@@ -435,7 +435,7 @@ def test_measurement_basis_list(backend):
     c.add(gates.H(2))
     c.add(gates.X(3))
     c.add(gates.M(0, 1, 2, 3, basis=[gates.X, gates.Z, gates.X, gates.Z]))
-    result = c(nshots=100)
+    result = backend.execute_circuit(c, nshots=100)
     assert result.frequencies() == {"0011": 100}
     assert (
         c.draw()

tests/test_backends.py

@@ -253,7 +253,7 @@ def test_collapse_error(backend):
     c = models.Circuit(1)
     m = c.add(gates.M(0, collapse=True))
     with pytest.raises(Exception) as exc_info:
-        c()
+        backend.execute_circuit(c)
     assert (
         str(exc_info.value)
         == "The circuit contains only collapsing measurements (`collapse=True`) but `density_matrix=False`. Please set `density_matrix=True` to retrieve the final state after execution."

tests/test_backends_qulacs.py

@@ -47,26 +47,32 @@ def test_circuit_init():
     assert c.nqubits == 2
 
 
-def test_eigenstate():
+def test_eigenstate(backend):
     nqubits = 3
     c = Circuit(nqubits)
     c.add(gates.M(*list(range(nqubits))))
     c2 = initialize(nqubits, eigenstate="-")
-    assert c(nshots=100, initial_state=c2).frequencies() == {"111": 100}
+    assert backend.execute_circuit(c, nshots=100, initial_state=c2).frequencies() == {
+        "111": 100
+    }
     c2 = initialize(nqubits, eigenstate="+")
-    assert c(nshots=100, initial_state=c2).frequencies() == {"000": 100}
+    assert backend.execute_circuit(c, nshots=100, initial_state=c2).frequencies() == {
+        "000": 100
+    }
 
     with pytest.raises(NotImplementedError):
         c2 = initialize(nqubits, eigenstate="x")
 
 
-def test_initialize():
+def test_initialize(backend):
     nqubits = 3
     for gate in [gates.X, gates.Y, gates.Z]:
         c = Circuit(nqubits)
         c.add(gates.M(*list(range(nqubits)), basis=gate))
         c2 = initialize(nqubits, basis=gate)
-        assert c(nshots=100, initial_state=c2).frequencies() == {"000": 100}
+        assert backend.execute_circuit(
+            c, nshots=100, initial_state=c2
+        ).frequencies() == {"000": 100}
 
 
 @pytest.mark.parametrize("nqubits", [0, -10, 2.5])

tests/test_measurements.py

@@ -7,7 +7,7 @@
 from qibo.models.iqae import IQAE
 
 
-def test_iqae_init(backend=None):
+def test_iqae_init(backend):
     A = Circuit(3 + 1)
     Q = Circuit(3 + 1)
     alpha = 0.05
@@ -23,7 +23,7 @@ def test_iqae_init(backend=None):
     assert iqae.method == method
 
 
-def test_iqae_init_raising_errors(backend=None):
+def test_iqae_init_raising_errors(backend):
     A = Circuit(3 + 1)
     Q = Circuit(4 + 1)
     # incorrect values of:
@@ -50,7 +50,7 @@ def test_iqae_init_raising_errors(backend=None):
     results = iqae.execute(backend=backend)
 
 
-def test_iqae_execution(backend=None):
+def test_iqae_execution(backend):
     # Let's check if we get the correct result for the integral of Sin(x)^2 from 0 to 1
     nbit = 3
     A = A_circ(qx=list(range(nbit)), qx_measure=nbit, nbit=nbit, b_max=1, b_min=0)

tests/test_measurements_collapse.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 import qibo
-from qibo import gates
+from qibo import gates, set_backend
 from qibo.models import Circuit
 from qibo.models.qcnn import QuantumCNN
 
@@ -13,6 +13,7 @@
 
 def test_classifier_circuit2():
     """ """
+    set_backend("numpy")
     nqubits = 2
     nlayers = int(nqubits / 2)
     init_state = np.ones(2**nqubits) / np.sqrt(2**nqubits)  #
@@ -77,6 +78,7 @@ def get_real_vector2():
 
 def test_classifier_circuit4():
     """ """
+    set_backend("numpy")
     nqubits = 4
     nlayers = int(nqubits / 2)
     init_state = np.ones(2**nqubits) / np.sqrt(2**nqubits)  #
@@ -278,6 +280,8 @@ def test_1_qubit_classifier_circuit_error():
 def test_qcnn_training():
     import random
 
+    set_backend("numpy")
+
     # generate 2 random states and labels for pytest
     data = np.zeros([2, 16])
     for i in range(2):