superoperator transformations

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

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

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

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)
+