From 3889f115180e2dad00f616d2405bc79e43a510fa Mon Sep 17 00:00:00 2001 From: Alejandro Sopena Date: Mon, 15 Apr 2024 18:46:27 +0200 Subject: [PATCH] superoperator transformations --- src/qibojit/backends/cpu.py | 5 + src/qibojit/backends/gpu.py | 18 ++++ .../backends/utils_quantum_info_cpu.py | 96 +++++++++++++++++++ test_numba.py | 52 ++++++++++ 4 files changed, 171 insertions(+) create mode 100644 src/qibojit/backends/utils_quantum_info_cpu.py create mode 100644 test_numba.py diff --git a/src/qibojit/backends/cpu.py b/src/qibojit/backends/cpu.py index 4ae8deaf..995b00cd 100644 --- a/src/qibojit/backends/cpu.py +++ b/src/qibojit/backends/cpu.py @@ -303,3 +303,8 @@ def sample_frequencies(self, probabilities, nshots): # def calculate_frequencies(self, samples): Inherited from ``NumpyBackend`` # def assert_allclose(self, value, target, rtol=1e-7, atol=0.0): Inherited from ``NumpyBackend`` + + def _load_quantum_info_utils(self): + from qibojit.backends import utils_quantum_info_cpu + + return utils_quantum_info_cpu \ No newline at end of file diff --git a/src/qibojit/backends/gpu.py b/src/qibojit/backends/gpu.py index c00d233c..da455772 100644 --- a/src/qibojit/backends/gpu.py +++ b/src/qibojit/backends/gpu.py @@ -509,6 +509,12 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): expd = self.cp.diag(self.cp.exp(-1j * a * eigenvalues)) ud = self.cp.transpose(self.cp.conj(eigenvectors)) return self.cp.matmul(eigenvectors, self.cp.matmul(expd, ud)) + + @staticmethod + def _load_quantum_info_utils(): + from qibo.quantum_info import _utils_numpy + + return _utils_numpy class CuQuantumBackend(CupyBackend): # pragma: no cover @@ -804,6 +810,12 @@ def collapse_state(self, state, qubits, shot, nqubits, normalize=True): state = state / norm return state + + @staticmethod + def _load_quantum_info_utils(): + from qibo.quantum_info import _utils_numpy + + return _utils_numpy class MultiGpuOps: # pragma: no cover @@ -901,3 +913,9 @@ def revert_swaps(self, pieces, swap_pairs): q1, q2 = q2, q1 pieces = self.swap(pieces, q1, q2) return pieces + + @staticmethod + def _load_quantum_info_utils(): + from qibo.quantum_info import _utils_numpy + + return _utils_numpy diff --git a/src/qibojit/backends/utils_quantum_info_cpu.py b/src/qibojit/backends/utils_quantum_info_cpu.py new file mode 100644 index 00000000..b4f29bc9 --- /dev/null +++ b/src/qibojit/backends/utils_quantum_info_cpu.py @@ -0,0 +1,96 @@ +import numpy as np +from numba import njit + + +@njit(cache=True) +def _dim_shape_vectorization(state): + dim = len(state) + nqubits = int(np.log2(dim)) + shape = [2] * 2 * nqubits + new_axis = [i for qubit in range(nqubits) for i in (qubit + nqubits, qubit)] + return shape, new_axis + + +@njit(parallel=True,cache=True) +def vectorization_row_colum(state, order): + if order == "row": + state1 = np.ravel(state) + return state1 + elif order == "column": + state1 = np.transpose(np.ravel(state)) + return state1 + else: + raise ValueError("Invalid order. Use 'row' or 'column'.") + + +@njit(parallel=True,cache=True) +def vectorization_system(state, shape=None, axis=None): + state1 = np.reshape(state, shape) + state1 = np.transpose(state1, axes=axis) + state1 = np.ravel(state1) + return state1 + +@njit(parallel=True,cache=True) +def outer_conj(state): + return np.outer(state, np.conj(state)) + +def vectorization(state, order): + if len(state.shape) == 1: + state = outer_conj(state) + if order == "system": + shape, new_axis = _dim_shape_vectorization(state) + shape = tuple(shape) + new_axis = tuple(new_axis) + state = vectorization_system(state, shape, new_axis) + else: + state = vectorization_row_colum(state, order) + return state + +@njit(parallel=True,cache=True) +def unvectorization_system(state, shape1, shape2, axis): + state1 = np.reshape(state, shape1) + state1 = np.copy(np.transpose(state1, axes=axis)) + state1 = np.reshape(state1, shape2) + return state1 + +def unvectorization(state, order): + dim = int(np.sqrt(len(state))) + if order in ["row", "column"]: + order = "C" if order == "row" else "F" + state1 = np.reshape(state, (dim, dim), order=order) + return state1 + else: + state1 = state + nqubits = int(np.log2(dim)) + shape1 = tuple([2] * 2 * nqubits) + shape2 = tuple([2**nqubits] * 2) + axis = list(np.arange(0, 2 * nqubits)) + axis = tuple(axis[1::2] + axis[0::2]) + + state1 = unvectorization_system(state1, shape1, shape2, axis) + return state1 + + +def to_choi(channel, order): + channel = vectorization(channel, order=order) + channel = outer_conj(channel) #np.outer(channel, np.conj(channel)) + + return channel + +def to_liouville(channel, order): + channel = to_choi(channel, order=order) + channel = _reshuffling(channel, order=order) + + return channel + +def _reshuffling(super_op, order): + dim = np.sqrt(super_op.shape[0]) + dim = int(dim) + super_op = np.reshape(super_op, [dim] * 4) + + axes = [1, 2] if order == "row" else [0, 3] + super_op = np.swapaxes(super_op, *axes) + + super_op = np.reshape(super_op, [dim**2, dim**2]) + + return super_op \ No newline at end of file diff --git a/test_numba.py b/test_numba.py new file mode 100644 index 00000000..7c4e13ac --- /dev/null +++ b/test_numba.py @@ -0,0 +1,52 @@ +# from qibojit.backends.utils_quantum_info_cpu import vectorization +import numpy as np +from qibo.backends import set_backend, construct_backend +set_backend("numpy") + + + + +# vectorization(state, order) + +# print(vectorization.signatures) + +from time import time +from qibo.quantum_info.superoperator_transformations import vectorization, unvectorization +# from qibojit.backends.utils_quantum_info_cpu import vectorization, get_dim_shape + +# backend = construct_backend("qibojit",platform="numba") + + +state = np.random.rand(2**28) +order = "system" + + +backend = construct_backend("numpy") +utils = backend._load_quantum_info_utils() +utils.unvectorization +t = time() +unvectorization(state, order, backend=backend) +print("Elapsed time:", time() - t) +t = time() +unvectorization(state, order, backend=backend) +print("Elapsed time:", time() - t) + +backend = construct_backend("qibojit",platform="numba") +t = time() +unvectorization(state, order, backend=backend) +print("Elapsed time:", time() - t) + +t = time() +unvectorization(state, order, backend=backend) +print("Elapsed time:", time() - t) +# dim = len(state) +# nqubits = int(np.log2(dim)) +# shape = tuple([2] * 2 * nqubits) + +# new_axis = [int(x) for x in range(0)] +# for qubit in range(nqubits): +# new_axis += [qubit + nqubits, qubit] +# new_axis = tuple(new_axis) + +# vectorization(state, order,shape,new_axis) +