Cupy matrices

pyproject.toml

@@ -56,7 +56,7 @@ 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"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git", branch="cupy_matrices" }
 
 [tool.poetry.group.tests]
 optional = true
@@ -70,7 +70,7 @@ pytest-cov = "^4.0.0"
 pylint = "3.1.0"
 matplotlib = "^3.7.0"
 torch = "^2.1.1,<2.4"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git", branch="cupy_matrices" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 qiboml = { git  = "https://github.com/qiboteam/qiboml.git" }
 stim = "^1.12.0"
@@ -89,19 +89,19 @@ test-docs = "make -C doc doctest"
 optional = true
 
 [tool.poetry.group.cuda11.dependencies]
-cupy-cuda11x = "^12.0.0"
+cupy-cuda11x = "^13.1.0"
 cuquantum-python-cu11 = "^23.3.0"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git", branch="cupy_matrices" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 
 
 [tool.poetry.group.cuda12]
 optional = true
 
 [tool.poetry.group.cuda12.dependencies]
-cupy-cuda12x = "^12.0.0"
+cupy-cuda12x = "^13.1.0"
 cuquantum-python-cu12 = "^23.3.0"
-qibojit = { git = "https://github.com/qiboteam/qibojit.git" }
+qibojit = { git = "https://github.com/qiboteam/qibojit.git", branch="cupy_matrices" }
 qibotn = { git = "https://github.com/qiboteam/qibotn.git" }
 
 [tool.poetry.extras]

src/qibo/backends/clifford.py

@@ -48,6 +48,10 @@ def __init__(self, engine=None):
         if engine == "numpy":
             pass
         elif engine == "numba":
+            from numba import set_num_threads
+
+            set_num_threads(1)
+
             from qibojit.backends import (  # pylint: disable=C0415
                 clifford_operations_cpu,
             )

src/qibo/backends/npmatrices.py

@@ -15,7 +15,9 @@ def __init__(self, dtype):
         self.np = np
 
     def _cast(self, x, dtype):
-        return self.np.array(x, dtype=dtype)
+        if isinstance(x, list):
+            return self.np.array(x, dtype=dtype)
+        return x.astype(dtype)
 
     @cached_property
     def H(self):

src/qibo/optimizers.py

@@ -84,7 +84,9 @@ def myloss(parameters, circuit):
                 RuntimeError,
                 "The keyword 'bounds' cannot be used with the cma optimizer. Please use 'options' instead as defined by the cma documentation: ex. options['bounds'] = [0.0, 1.0].",
             )
-        return cmaes(loss, initial_parameters, args, callback, options)
+        return cmaes(
+            loss, backend.to_numpy(initial_parameters), args, callback, options
+        )
     elif method == "sgd":
         from qibo.backends import _check_backend
 
@@ -208,13 +210,13 @@ def newtonian(
             callback=callback,
             options=options,
         )
-        m = o.run(initial_parameters)
+        m = o.run(backend.to_numpy(initial_parameters))
     else:
         from scipy.optimize import minimize
 
         m = minimize(
             loss,
-            initial_parameters,
+            backend.to_numpy(initial_parameters),
             args=args,
             method=method,
             jac=jac,

src/qibo/tomography/gate_set_tomography.py

@@ -187,7 +187,7 @@ def _gate_tomography(
         gate = gate.__class__(*gate.qubits, **gate.init_kwargs)
 
     # GST for empty circuit or with gates
-    matrix_jk = np.zeros((4**nqubits, 4**nqubits))
+    matrix_jk = 1j * np.zeros((4**nqubits, 4**nqubits))
     for k in range(4**nqubits):
         circ = Circuit(nqubits, density_matrix=True)
         circ.add(_prepare_state(k, nqubits))
@@ -211,7 +211,7 @@ def _gate_tomography(
                     backend.execute_circuit(new_circ, nshots=nshots).frequencies()
                 )
             matrix_jk[j, k] = exp_val
-    return matrix_jk
+    return backend.cast(matrix_jk, dtype=matrix_jk.dtype)
 
 
 def GST(
@@ -327,20 +327,20 @@ def GST(
             if np.linalg.det(gauge_matrix) == 0:
                 raise_error(ValueError, "Matrix is not invertible")
         else:
-            gauge_matrix = np.array(
+            gauge_matrix = backend.cast(
                 [[1, 1, 1, 1], [0, 0, 1, 0], [0, 0, 0, 1], [1, -1, 0, 0]]
             )
         PL_matrices = []
         gauge_matrix_1q = gauge_matrix
-        gauge_matrix_2q = np.kron(gauge_matrix, gauge_matrix)
+        gauge_matrix_2q = backend.np.kron(gauge_matrix, gauge_matrix)
         for matrix in matrices:
             gauge_matrix = gauge_matrix_1q if matrix.shape[0] == 4 else gauge_matrix_2q
             empty = empty_matrices[0] if matrix.shape[0] == 4 else empty_matrices[1]
             PL_matrices.append(
                 gauge_matrix
-                @ np.linalg.inv(empty)
+                @ backend.np.linalg.inv(empty)
                 @ matrix
-                @ np.linalg.inv(gauge_matrix)
+                @ backend.np.linalg.inv(gauge_matrix)
             )
         matrices = PL_matrices
 

src/qibo/transpiler/unitary_decompositions.py

@@ -32,8 +32,8 @@ def u3_decomposition(unitary, backend):
     minus = backend.np.angle(su2[1, 0])
     phi = plus + minus
     lam = plus - minus
-
-    return theta, phi, lam
+    # explicit conversion to float to avoid issue on GPU
+    return float(theta), float(phi), float(lam)
 
 
 def calculate_psi(unitary, backend, magic_basis=magic_basis):
@@ -51,16 +51,6 @@ def calculate_psi(unitary, backend, magic_basis=magic_basis):
     Returns:
         ndarray: Eigenvectors in the computational basis and eigenvalues of :math:`U^{T} U`.
     """
-
-    if backend.__class__.__name__ in [
-        "CupyBackend",
-        "CuQuantumBackend",
-    ]:  # pragma: no cover
-        raise_error(
-            NotImplementedError,
-            f"{backend.__class__.__name__} does not support `linalg.eig.`",
-        )
-
     magic_basis = backend.cast(magic_basis)
     unitary = backend.cast(unitary)
     # write unitary in magic basis
@@ -71,10 +61,12 @@ def calculate_psi(unitary, backend, magic_basis=magic_basis):
     )
     # construct and diagonalize UT_U
     ut_u = backend.np.transpose(u_magic, (1, 0)) @ u_magic
-    if backend.__class__.__name__ != "PyTorchBackend":
+    if backend.__class__.__name__ not in ("PyTorchBackend", "TensorflowBackend"):
         # eig seems to have a different behavior based on backend/hardware,
         # use np.round to increase precision seems to fix the issue
-        eigvals, psi_magic = np.linalg.eig(np.round(ut_u, decimals=20))
+        eigvals, psi_magic = backend.calculate_eigenvectors(
+            np.round(ut_u, decimals=20), hermitian=False
+        )
     else:
         eigvals, psi_magic = backend.calculate_eigenvectors(ut_u, hermitian=False)
     # orthogonalize eigenvectors in the case of degeneracy (Gram-Schmidt)

tests/test_backends.py

@@ -7,6 +7,8 @@
 from qibo import construct_backend, gates, list_available_backends, set_backend
 from qibo.backends import MetaBackend
 
+from .conftest import AVAILABLE_BACKENDS
+
 ####################### Test `matrix` #######################
 GATES = [
     ("H", (0,), np.array([[1, 1], [1, -1]]) / np.sqrt(2)),
@@ -119,7 +121,10 @@ def test_list_available_backends():
     available_backends = {
         "numpy": True,
         "qulacs": qulacs,
-        "qibojit": {"numba": True, "cupy": False, "cuquantum": False},
+        "qibojit": {
+            platform: any(platform in backend for backend in AVAILABLE_BACKENDS)
+            for platform in ["numba", "cupy", "cuquantum"]
+        },
         "qibolab": False,
         "qibo-cloud-backends": False,
         "qibotn": {"cutensornet": False, "qutensornet": True},

tests/test_backends_clifford.py

@@ -4,7 +4,6 @@
 
 import numpy as np
 import pytest
-from qiboml.backends import PyTorchBackend, TensorflowBackend
 
 from qibo import Circuit, gates, get_backend, set_backend
 from qibo.backends import CliffordBackend, NumpyBackend
@@ -16,9 +15,10 @@
 
 
 def construct_clifford_backend(backend):
-    if (
-        isinstance(backend, (TensorflowBackend, PyTorchBackend))
-        or backend.__class__.__name__ == "CuQuantumBackend"
+    if backend.__class__.__name__ in (
+        "TensorflowBackend",
+        "PyTorchBackend",
+        "CuQuantumBackend",
     ):
         with pytest.raises(NotImplementedError):
             clifford_backend = CliffordBackend(backend.name)

tests/test_backends_qibotn.py

@@ -1,11 +1,7 @@
-import os
-
-from qibo import get_backend, get_threads, set_backend
-
-# Force quimb to use qibojit default number of threads.
-os.environ["NUMBA_NUM_THREADS"] = f"{get_threads()}"
 from qibotn.backends.quimb import QuimbBackend
 
+from qibo import get_backend, set_backend
+
 
 def test_backend_qibotn():
     set_backend(backend="qibotn", platform="qutensornet", runcard=None)

tests/test_models_variational.py

@@ -144,7 +144,8 @@ def test_vqe(backend, method, options, compile, filename):
 
     def callback(parameters, loss_values=loss_values, vqe=v):
         vqe.circuit.set_parameters(parameters)
-        loss_values.append(vqe.hamiltonian.expectation(vqe.circuit().state()))
+        state = vqe.backend.execute_circuit(vqe.circuit).state()
+        loss_values.append(vqe.hamiltonian.expectation(state))
 
     best, params, _ = v.minimize(
         initial_parameters,

tests/test_quantum_info_clifford.py

@@ -3,7 +3,6 @@
 
 import numpy as np
 import pytest
-from qiboml.backends import PyTorchBackend, TensorflowBackend
 
 from qibo import Circuit, gates, matrices
 from qibo.backends import CliffordBackend
@@ -18,9 +17,10 @@
 
 
 def construct_clifford_backend(backend):
-    if (
-        isinstance(backend, (TensorflowBackend, PyTorchBackend))
-        or backend.__class__.__name__ == "CuQuantumBackend"
+    if backend.__class__.__name__ in (
+        "TensorflowBackend",
+        "PyTorchBackend",
+        "CuQuantumBackend",
     ):
         with pytest.raises(NotImplementedError):
             clifford_backend = CliffordBackend(backend.name)

tests/test_quantum_info_entropies.py

@@ -436,17 +436,8 @@ def test_von_neumann_entropy(backend, base, check_hermitian):
             state, base=base, check_hermitian="False", backend=backend
         )
 
-    if backend.__class__.__name__ in ["CupyBackend", "CuQuantumBackend"]:
-        with pytest.raises(NotImplementedError):
-            state = random_unitary(4, backend=backend)
-            test = von_neumann_entropy(
-                state, base=base, check_hermitian=True, backend=backend
-            )
-    else:
-        state = random_unitary(4, backend=backend)
-        test = von_neumann_entropy(
-            state, base=base, check_hermitian=True, backend=backend
-        )
+    state = random_unitary(4, backend=backend)
+    test = von_neumann_entropy(state, base=base, check_hermitian=True, backend=backend)
 
     state = np.array([1.0, 0.0])
     state = backend.cast(state, dtype=state.dtype)
@@ -842,16 +833,6 @@ def test_entanglement_entropy(backend, bipartition, base, check_hermitian):
             check_hermitian=check_hermitian,
             backend=backend,
         )
-    if backend.__class__.__name__ == "CupyBackend":
-        with pytest.raises(NotImplementedError):
-            state = random_unitary(4, backend=backend)
-            test = entanglement_entropy(
-                state,
-                bipartition=bipartition,
-                base=base,
-                check_hermitian=True,
-                backend=backend,
-            )
 
     # Bell state
     state = np.array([1.0, 0.0, 0.0, 1.0]) / np.sqrt(2)

tests/test_tomography_gate_set_tomography.py

@@ -216,8 +216,10 @@ def test_gate_tomography_noise_model(backend):
 )
 @pytest.mark.parametrize("pauli_liouville", [False, True])
 def test_GST(backend, target_gates, pauli_liouville):
-    T = np.array([[1, 1, 1, 1], [0, 0, 1, 0], [0, 0, 0, 1], [1, -1, 0, 0]])
-    T = backend.cast(T, dtype=T.dtype)
+    T = np.array(
+        [[1.0, 1, 1, 1], [0, 0, 1, 0], [0, 0, 0, 1], [1, -1, 0, 0]], dtype=np.complex128
+    )
+    T = backend.cast(T)
     target_matrices = [g.matrix(backend=backend) for g in target_gates]
     # superoperator representation of the target gates in the Pauli basis
     target_matrices = [
@@ -237,11 +239,11 @@ def test_GST(backend, target_gates, pauli_liouville):
             pauli_liouville=pauli_liouville,
             backend=backend,
         )
-        T_2q = np.kron(T, T)
+        T_2q = backend.np.kron(T, T)
         for target, estimate in zip(target_matrices, approx_gates):
             if not pauli_liouville:
                 G = empty_1q if estimate.shape[0] == 4 else empty_2q
-                G_inv = np.linalg.inv(G)
+                G_inv = backend.np.linalg.inv(G)
                 T_matrix = T if estimate.shape[0] == 4 else T_2q
                 estimate = T_matrix @ G_inv @ estimate @ G_inv
             backend.assert_allclose(