diff --git a/src/qibo/backends/npmatrices.py b/src/qibo/backends/npmatrices.py index fb7b274906..58e78f4ca1 100644 --- a/src/qibo/backends/npmatrices.py +++ b/src/qibo/backends/npmatrices.py @@ -1,3 +1,5 @@ +import cmath +import math from functools import cached_property from qibo.config import raise_error @@ -15,9 +17,13 @@ def __init__(self, dtype): def _cast(self, x, dtype): return self.np.array(x, dtype=dtype) + # This method is used to cast the parameters of the gates to the right type for other backends + def _cast_parameter(self, x): + return x + @cached_property def H(self): - return self._cast([[1, 1], [1, -1]], dtype=self.dtype) / self.np.sqrt(2) + return self._cast([[1, 1], [1, -1]], dtype=self.dtype) / math.sqrt(2) @cached_property def X(self): @@ -50,19 +56,17 @@ def SDG(self): @cached_property def T(self): return self._cast( - [[1, 0], [0, self.np.exp(1j * self.np.pi / 4.0)]], dtype=self.dtype + [[1 + 0j, 0], [0, cmath.exp(1j * math.pi / 4.0)]], dtype=self.dtype ) @cached_property def TDG(self): return self._cast( - [[1, 0], [0, self.np.exp(-1j * self.np.pi / 4.0)]], dtype=self.dtype + [[1 + 0j, 0], [0, cmath.exp(-1j * math.pi / 4.0)]], dtype=self.dtype ) def I(self, n=2): - # dtype=complex is necessary for pytorch backend, - # _cast will take care of casting in the right dtype for all the backends - return self._cast(self.np.eye(n, dtype=complex), dtype=self.dtype) + return self._cast(self.np.eye(n), dtype=self.dtype) def Align(self, delay, n=2): return self._cast(self.I(n), dtype=self.dtype) @@ -71,20 +75,25 @@ def M(self): # pragma: no cover raise_error(NotImplementedError) def RX(self, theta): + theta = self._cast_parameter(theta) cos = self.np.cos(theta / 2.0) + 0j isin = -1j * self.np.sin(theta / 2.0) return self._cast([[cos, isin], [isin, cos]], dtype=self.dtype) def RY(self, theta): + theta = self._cast_parameter(theta) cos = self.np.cos(theta / 2.0) + 0j sin = self.np.sin(theta / 2.0) + 0j return self._cast([[cos, -sin], [sin, cos]], dtype=self.dtype) def RZ(self, theta): + theta = self._cast_parameter(theta) phase = self.np.exp(0.5j * theta) return self._cast([[self.np.conj(phase), 0], [0, phase]], dtype=self.dtype) def PRX(self, theta, phi): + theta = self._cast_parameter(theta) + phi = self._cast_parameter(phi) cos = self.np.cos(theta / 2) sin = self.np.sin(theta / 2) exponent1 = -1.0j * self.np.exp(-1.0j * phi) @@ -95,29 +104,36 @@ def PRX(self, theta, phi): ) def GPI(self, phi): + phi = self._cast_parameter(phi) phase = self.np.exp(1.0j * phi) return self._cast([[0, self.np.conj(phase)], [phase, 0]], dtype=self.dtype) def GPI2(self, phi): + phi = self._cast_parameter(phi) phase = self.np.exp(1.0j * phi) return self._cast( [[1, -1.0j * self.np.conj(phase)], [-1.0j * phase, 1]], dtype=self.dtype - ) / self.np.sqrt(2) + ) / math.sqrt(2) def U1(self, theta): + theta = self._cast_parameter(theta) phase = self.np.exp(1j * theta) return self._cast([[1, 0], [0, phase]], dtype=self.dtype) def U2(self, phi, lam): + phi = self._cast_parameter(phi) + lam = self._cast_parameter(lam) eplus = self.np.exp(1j * (phi + lam) / 2.0) eminus = self.np.exp(1j * (phi - lam) / 2.0) return self._cast( - [[self.np.conj(eplus), -self.np.conj(eminus)], [eminus, eplus]] - / self.np.sqrt(2), + [[self.np.conj(eplus), -self.np.conj(eminus)], [eminus, eplus]], dtype=self.dtype, - ) + ) / math.sqrt(2) def U3(self, theta, phi, lam): + theta = self._cast_parameter(theta) + phi = self._cast_parameter(phi) + lam = self._cast_parameter(lam) cost = self.np.cos(theta / 2) sint = self.np.sin(theta / 2) eplus = self.np.exp(1j * (phi + lam) / 2.0) @@ -131,8 +147,10 @@ def U3(self, theta, phi, lam): ) def U1q(self, theta, phi): + theta = self._cast_parameter(theta) + phi = self._cast_parameter(phi) return self._cast( - self.U3(theta, phi - self.np.pi / 2, self.np.pi / 2 - phi), dtype=self.dtype + self.U3(theta, phi - math.pi / 2, math.pi / 2 - phi), dtype=self.dtype ) @cached_property @@ -161,7 +179,7 @@ def CZ(self): @cached_property def CSX(self): - a = (1 + 1j) / 2 + a = self._cast_parameter((1 + 1j) / 2) b = self.np.conj(a) return self._cast( [ @@ -175,7 +193,7 @@ def CSX(self): @cached_property def CSXDG(self): - a = (1 - 1j) / 2 + a = self._cast_parameter((1 - 1j) / 2) b = self.np.conj(a) return self._cast( [ @@ -188,6 +206,7 @@ def CSXDG(self): ) def CRX(self, theta): + theta = self._cast_parameter(theta) cos = self.np.cos(theta / 2.0) + 0j isin = -1j * self.np.sin(theta / 2.0) matrix = [ @@ -199,12 +218,14 @@ def CRX(self, theta): return self._cast(matrix, dtype=self.dtype) def CRY(self, theta): + theta = self._cast_parameter(theta) cos = self.np.cos(theta / 2.0) + 0j sin = self.np.sin(theta / 2.0) + 0j matrix = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, cos, -sin], [0, 0, sin, cos]] return self._cast(matrix, dtype=self.dtype) def CRZ(self, theta): + theta = self._cast_parameter(theta) phase = self.np.exp(0.5j * theta) matrix = [ [1, 0, 0, 0], @@ -215,6 +236,7 @@ def CRZ(self, theta): return self._cast(matrix, dtype=self.dtype) def CU1(self, theta): + theta = self._cast_parameter(theta) phase = self.np.exp(1j * theta) matrix = [ [1, 0, 0, 0], @@ -225,8 +247,10 @@ def CU1(self, theta): return self._cast(matrix, dtype=self.dtype) def CU2(self, phi, lam): - eplus = self.np.exp(1j * (phi + lam) / 2.0) / self.np.sqrt(2) - eminus = self.np.exp(1j * (phi - lam) / 2.0) / self.np.sqrt(2) + phi = self._cast_parameter(phi) + lam = self._cast_parameter(lam) + eplus = self.np.exp(1j * (phi + lam) / 2.0) / math.sqrt(2) + eminus = self.np.exp(1j * (phi - lam) / 2.0) / math.sqrt(2) matrix = [ [1, 0, 0, 0], [0, 1, 0, 0], @@ -236,6 +260,9 @@ def CU2(self, phi, lam): return self._cast(matrix, dtype=self.dtype) def CU3(self, theta, phi, lam): + theta = self._cast_parameter(theta) + phi = self._cast_parameter(phi) + lam = self._cast_parameter(lam) cost = self.np.cos(theta / 2) sint = self.np.sin(theta / 2) eplus = self.np.exp(1j * (phi + lam) / 2.0) @@ -271,8 +298,8 @@ def SiSWAP(self): return self._cast( [ [1 + 0j, 0j, 0j, 0j], - [0j, 1 / self.np.sqrt(2) + 0j, 1j / self.np.sqrt(2), 0j], - [0j, 1j / self.np.sqrt(2), 1 / self.np.sqrt(2) + 0j, 0j], + [0j, 1 / math.sqrt(2) + 0j, 1j / math.sqrt(2), 0j], + [0j, 1j / math.sqrt(2), 1 / math.sqrt(2) + 0j, 0j], [0j, 0j, 0j, 1 + 0j], ], dtype=self.dtype, @@ -283,8 +310,8 @@ def SiSWAPDG(self): return self._cast( [ [1 + 0j, 0j, 0j, 0j], - [0j, 1 / self.np.sqrt(2) + 0j, -1j / self.np.sqrt(2), 0j], - [0j, -1j / self.np.sqrt(2), 1 / self.np.sqrt(2) + 0j, 0j], + [0j, 1 / math.sqrt(2) + 0j, -1j / math.sqrt(2), 0j], + [0j, -1j / math.sqrt(2), 1 / math.sqrt(2) + 0j, 0j], [0j, 0j, 0j, 1 + 0j], ], dtype=self.dtype, @@ -297,6 +324,8 @@ def FSWAP(self): ) def fSim(self, theta, phi): + theta = self._cast_parameter(theta) + phi = self._cast_parameter(phi) cost = self.np.cos(theta) + 0j isint = -1j * self.np.sin(theta) phase = self.np.exp(-1j * phi) @@ -312,12 +341,12 @@ def fSim(self, theta, phi): @cached_property def SYC(self): - cost = self.np.cos(self.np.pi / 2) + 0j - isint = -1j * self.np.sin(self.np.pi / 2) - phase = self.np.exp(-1j * self.np.pi / 6) + cost = math.cos(math.pi / 2) + 0j + isint = -1j * math.sin(math.pi / 2) + phase = cmath.exp(-1j * math.pi / 6) return self._cast( [ - [1, 0, 0, 0], + [1 + 0j, 0, 0, 0], [0, cost, isint, 0], [0, isint, cost, 0], [0, 0, 0, phase], @@ -326,6 +355,7 @@ def SYC(self): ) def GeneralizedfSim(self, u, phi): + phi = self._cast_parameter(phi) phase = self.np.exp(-1j * phi) return self._cast( [ @@ -338,6 +368,7 @@ def GeneralizedfSim(self, u, phi): ) def RXX(self, theta): + theta = self._cast_parameter(theta) cos = self.np.cos(theta / 2.0) + 0j isin = -1j * self.np.sin(theta / 2.0) return self._cast( @@ -351,6 +382,7 @@ def RXX(self, theta): ) def RYY(self, theta): + theta = self._cast_parameter(theta) cos = self.np.cos(theta / 2.0) + 0j isin = -1j * self.np.sin(theta / 2.0) return self._cast( @@ -364,6 +396,7 @@ def RYY(self, theta): ) def RZZ(self, theta): + theta = self._cast_parameter(theta) phase = self.np.exp(0.5j * theta) return self._cast( [ @@ -376,6 +409,7 @@ def RZZ(self, theta): ) def RZX(self, theta): + theta = self._cast_parameter(theta) cos, sin = self.np.cos(theta / 2) + 0j, self.np.sin(theta / 2) + 0j return self._cast( [ @@ -388,6 +422,7 @@ def RZX(self, theta): ) def RXXYY(self, theta): + theta = self._cast_parameter(theta) cos, sin = self.np.cos(theta / 2) + 0j, self.np.sin(theta / 2) + 0j return self._cast( [ @@ -400,6 +435,11 @@ def RXXYY(self, theta): ) def MS(self, phi0, phi1, theta): + phi0, phi1, theta = ( + self._cast_parameter(phi0), + self._cast_parameter(phi1), + self._cast_parameter(theta), + ) plus = self.np.exp(1.0j * (phi0 + phi1)) minus = self.np.exp(1.0j * (phi0 - phi1)) cos = self.np.cos(theta / 2) + 0j @@ -415,6 +455,7 @@ def MS(self, phi0, phi1, theta): ) def GIVENS(self, theta): + theta = self._cast_parameter(theta) return self._cast( [ [1, 0, 0, 0], @@ -438,7 +479,7 @@ def ECR(self): [-1j, 1 + 0j, 0j, 0j], ], dtype=self.dtype, - ) / self.np.sqrt(2) + ) / math.sqrt(2) @cached_property def TOFFOLI(self): @@ -473,6 +514,7 @@ def CCZ(self): ) def DEUTSCH(self, theta): + theta = self._cast_parameter(theta) sin = self.np.sin(theta) + 0j # 0j necessary for right tensorflow dtype cos = self.np.cos(theta) + 0j return self._cast( diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index 5760a62bda..ec4329ff00 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -156,7 +156,6 @@ def matrix_fused(self, fgate): return self.cast(matrix.toarray()) def apply_gate(self, gate, state, nqubits): - state = self.cast(state) state = self.np.reshape(state, nqubits * (2,)) matrix = gate.matrix(self) if gate.is_controlled_by: @@ -718,23 +717,27 @@ def calculate_norm_density_matrix(self, state, order="nuc"): return self.np.linalg.norm(state, ord=order) def calculate_overlap(self, state1, state2): - return self.np.abs(self.np.sum(np.conj(self.cast(state1)) * self.cast(state2))) + return self.np.abs( + self.np.sum(self.np.conj(self.cast(state1)) * self.cast(state2)) + ) def calculate_overlap_density_matrix(self, state1, state2): return self.np.trace( self.np.matmul(self.np.conj(self.cast(state1)).T, self.cast(state2)) ) - def calculate_eigenvalues(self, matrix, k=6): + def calculate_eigenvalues(self, matrix, k=6, hermitian=True): if self.issparse(matrix): log.warning( "Calculating sparse matrix eigenvectors because " "sparse modules do not provide ``eigvals`` method." ) return self.calculate_eigenvectors(matrix, k=k)[0] - return np.linalg.eigvalsh(matrix) + if hermitian: + return np.linalg.eigvalsh(matrix) + return np.linalg.eigvals(matrix) - def calculate_eigenvectors(self, matrix, k=6): + def calculate_eigenvectors(self, matrix, k=6, hermitian=True): if self.issparse(matrix): if k < matrix.shape[0]: from scipy.sparse.linalg import eigsh @@ -742,7 +745,9 @@ def calculate_eigenvectors(self, matrix, k=6): return eigsh(matrix, k=k, which="SA") else: # pragma: no cover matrix = self.to_numpy(matrix) - return np.linalg.eigh(matrix) + if hermitian: + return np.linalg.eigh(matrix) + return np.linalg.eig(matrix) def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): if eigenvectors is None or self.issparse(matrix): diff --git a/src/qibo/backends/pytorch.py b/src/qibo/backends/pytorch.py index a8233d4e3b..69889f6f4e 100644 --- a/src/qibo/backends/pytorch.py +++ b/src/qibo/backends/pytorch.py @@ -8,17 +8,28 @@ class TorchMatrices(NumpyMatrices): - """Matrix representation of every gate as a torch Tensor.""" + """Matrix representation of every gate as a torch Tensor. - def __init__(self, dtype): + Args: + dtype (torch.dtype): Data type of the matrices. + requires_grad (bool): If ``True`` the matrices require gradient. + """ + + def __init__(self, dtype, requires_grad): import torch # pylint: disable=import-outside-toplevel super().__init__(dtype) - self.torch = torch + self.np = torch self.dtype = dtype + self.requires_grad = requires_grad def _cast(self, x, dtype): - return self.torch.as_tensor(x, dtype=dtype) + flattened = [item for sublist in x for item in sublist] + tensor_list = [self.np.as_tensor(i, dtype=dtype) for i in flattened] + return self.np.stack(tensor_list).reshape(len(x), len(x)) + + def _cast_parameter(self, x): + return self.np.tensor(x, dtype=self.dtype, requires_grad=self.requires_grad) def Unitary(self, u): return self._cast(u, dtype=self.dtype) @@ -29,6 +40,9 @@ def __init__(self): super().__init__() import torch # pylint: disable=import-outside-toplevel + # Global variable to enable or disable gradient calculation + self.gradients = True + self.np = torch self.name = "pytorch" @@ -39,16 +53,24 @@ def __init__(self): } self.dtype = self._torch_dtype(self.dtype) - self.matrices = TorchMatrices(self.dtype) + self.matrices = TorchMatrices(self.dtype, requires_grad=self.gradients) self.device = self.np.device("cuda:0" if torch.cuda.is_available() else "cpu") self.nthreads = 0 self.tensor_types = (self.np.Tensor, np.ndarray) # These functions in Torch works in a different way than numpy or have different names self.np.transpose = self.np.permute + self.np.copy = self.np.clone self.np.expand_dims = self.np.unsqueeze self.np.mod = self.np.remainder self.np.right_shift = self.np.bitwise_right_shift + self.np.sign = self.np.sgn + self.np.flatnonzero = lambda x: self.np.nonzero(x).flatten() + + def requires_grad(self, requires_grad): + """Enable or disable gradient calculation.""" + self.gradients = requires_grad + self.matrices.requires_grad = requires_grad def _torch_dtype(self, dtype): if dtype == "float": @@ -63,6 +85,7 @@ def cast( x, dtype=None, copy: bool = False, + requires_grad: bool = None, ): """Casts input as a Torch tensor of the specified dtype. @@ -77,7 +100,13 @@ def cast( Defaults to ``None``. copy (bool, optional): If ``True``, the input tensor is copied before casting. Defaults to ``False``. + requires_grad (bool, optional): If ``True``, the input tensor requires gradient. + If ``False``, the input tensor does not require gradient. + If ``None``, the default gradient setting of the backend is used. """ + if requires_grad is None: + requires_grad = self.gradients + if dtype is None: dtype = self.dtype elif isinstance(dtype, type): @@ -85,12 +114,16 @@ def cast( elif not isinstance(dtype, self.np.dtype): dtype = self._torch_dtype(str(dtype)) + # check if dtype is an integer to remove gradients + if dtype in [self.np.int32, self.np.int64, self.np.int8, self.np.int16]: + requires_grad = False + if isinstance(x, self.np.Tensor): x = x.to(dtype) elif isinstance(x, list) and all(isinstance(row, self.np.Tensor) for row in x): x = self.np.stack(x) else: - x = self.np.tensor(x, dtype=dtype) + x = self.np.tensor(x, dtype=dtype, requires_grad=requires_grad) if copy: return x.clone() @@ -143,11 +176,15 @@ def sample_shots(self, probabilities, nshots): self.cast(probabilities, dtype="float"), nshots, replacement=True ) - def calculate_eigenvalues(self, matrix, k=6): - return self.np.linalg.eigvalsh(matrix) # pylint: disable=not-callable + def calculate_eigenvalues(self, matrix, k=6, hermitian=True): + if hermitian: + return self.np.linalg.eigvalsh(matrix) # pylint: disable=not-callable + return self.np.linalg.eigvals(matrix) # pylint: disable=not-callable - def calculate_eigenvectors(self, matrix, k=6): - return self.np.linalg.eigh(matrix) # pylint: disable=not-callable + def calculate_eigenvectors(self, matrix, k=6, hermitian=True): + if hermitian: + return self.np.linalg.eigh(matrix) # pylint: disable=not-callable + return self.np.linalg.eig(matrix) # pylint: disable=not-callable def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): if eigenvectors is None or self.issparse(matrix): diff --git a/src/qibo/backends/tensorflow.py b/src/qibo/backends/tensorflow.py index 3d10da18f6..f45ec05ec9 100644 --- a/src/qibo/backends/tensorflow.py +++ b/src/qibo/backends/tensorflow.py @@ -19,6 +19,7 @@ def __init__(self, dtype): self.tf = tf self.np = tnp + self.np.linalg = tf.linalg def _cast(self, x, dtype): return self.tf.cast(x, dtype=dtype) @@ -42,6 +43,8 @@ def __init__(self): tnp.experimental_enable_numpy_behavior() self.tf = tf self.np = tnp + self.np.flatnonzero = np.flatnonzero + self.np.copy = np.copy self.versions = { "qibo": __version__, @@ -174,11 +177,15 @@ def calculate_norm_density_matrix(self, state, order="nuc"): return self.np.trace(state) return self.tf.norm(state, ord=order) - def calculate_eigenvalues(self, matrix, k=6): - return self.tf.linalg.eigvalsh(matrix) + def calculate_eigenvalues(self, matrix, k=6, hermitian=True): + if hermitian: + return self.tf.linalg.eigvalsh(matrix) + return self.tf.linalg.eigvals(matrix) - def calculate_eigenvectors(self, matrix, k=6): - return self.tf.linalg.eigh(matrix) + def calculate_eigenvectors(self, matrix, k=6, hermitian=True): + if hermitian: + return self.tf.linalg.eigh(matrix) + return self.tf.linalg.eig(matrix) def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): if eigenvectors is None or self.issparse(matrix): diff --git a/src/qibo/derivative.py b/src/qibo/derivative.py index 72772cc006..65d4ef87b0 100644 --- a/src/qibo/derivative.py +++ b/src/qibo/derivative.py @@ -1,6 +1,5 @@ import numpy as np -from qibo.backends.pytorch import PyTorchBackend from qibo.config import raise_error from qibo.hamiltonians.abstract import AbstractHamiltonian @@ -103,12 +102,6 @@ def circuit(nqubits = 1): # inheriting hamiltonian's backend backend = hamiltonian.backend - # TODO: make this work wih pytorch backend - if isinstance(backend, PyTorchBackend): - raise_error( - NotImplementedError, - "PyTorchBackend for the parameter shift rule is not supported.", - ) # getting the gate's type gate = circuit.associate_gates_with_parameters()[parameter_index] diff --git a/src/qibo/gates/channels.py b/src/qibo/gates/channels.py index f12df41a29..bde8c50273 100644 --- a/src/qibo/gates/channels.py +++ b/src/qibo/gates/channels.py @@ -97,7 +97,9 @@ def to_choi(self, nqubits: Optional[int] = None, order: str = "row", backend=Non kraus_op.append(gate) kraus_op = kraus_op.matrix(backend) kraus_op = vectorization(kraus_op, order=order, backend=backend) - super_op += coeff * np.outer(kraus_op, np.conj(kraus_op)) + super_op = super_op + coeff * backend.np.outer( + kraus_op, backend.np.conj(kraus_op) + ) del kraus_op return super_op @@ -176,7 +178,9 @@ def to_pauli_liouville( nqubits, normalize, pauli_order=pauli_order, backend=backend ) - super_op = unitary @ super_op @ np.transpose(np.conj(unitary)) + super_op = ( + unitary @ super_op @ backend.np.transpose(backend.np.conj(unitary), (1, 0)) + ) return super_op diff --git a/src/qibo/gates/gates.py b/src/qibo/gates/gates.py index 5b995b9e69..1dcc450d56 100644 --- a/src/qibo/gates/gates.py +++ b/src/qibo/gates/gates.py @@ -2508,24 +2508,12 @@ def __init__( } if check_unitary: - if unitary.__class__.__name__ == "Tensor": - import torch # pylint: disable=C0145 - - diag_function = torch.diag - all_function = torch.all - conj_function = torch.conj - transpose_function = torch.transpose - else: - diag_function = np.diag - all_function = np.all - conj_function = np.conj - transpose_function = np.transpose - - product = transpose_function(conj_function(unitary), (1, 0)) @ unitary - diagonals = all(np.abs(1 - diag_function(product)) < PRECISION_TOL) + engine = _check_engine(unitary) + product = engine.transpose(engine.conj(unitary), (1, 0)) @ unitary + diagonals = all(engine.abs(1 - engine.diag(product)) < PRECISION_TOL) off_diagonals = bool( - all_function( - np.abs(product - diag_function(diag_function(product))) + engine.all( + engine.abs(product - engine.diag(engine.diag(product))) < PRECISION_TOL ) ) @@ -2536,7 +2524,10 @@ def __init__( @Gate.parameters.setter def parameters(self, x): shape = self.parameters[0].shape - self._parameters = (np.reshape(x, shape),) + engine = _check_engine(self.parameters[0]) + # Reshape doesn't accept a tuple if engine is pytorch. + x = x[0] if type(x) is tuple else x + self._parameters = (engine.reshape(x, shape),) for gate in self.device_gates: # pragma: no cover gate.parameters = x @@ -2559,5 +2550,16 @@ def on_qubits(self, qubit_map): return gate def _dagger(self): - ud = np.conj(np.transpose(self.parameters[0])) + engine = _check_engine(self.parameters[0]) + ud = engine.conj(self.parameters[0].T) return self.__class__(ud, *self.target_qubits, **self.init_kwargs) + + +def _check_engine(array): + """Check if the array is a numpy or torch tensor and return the corresponding library.""" + if array.__class__.__name__ == "Tensor": + import torch # pylint: disable=C0415 + + return torch + else: + return np diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py index a2fc70d4ba..4bdec69d3b 100644 --- a/src/qibo/hamiltonians/hamiltonians.py +++ b/src/qibo/hamiltonians/hamiltonians.py @@ -6,7 +6,7 @@ import numpy as np import sympy -from qibo.backends import PyTorchBackend +from qibo.backends import PyTorchBackend, _check_backend from qibo.config import EINSUM_CHARS, log, raise_error from qibo.hamiltonians.abstract import AbstractHamiltonian from qibo.symbols import Z @@ -140,7 +140,12 @@ def expectation(self, state, normalize=False): def expectation_from_samples(self, freq, qubit_map=None): obs = self.matrix - if np.count_nonzero(obs - np.diag(np.diagonal(obs))) != 0: + if ( + self.backend.np.count_nonzero( + obs - self.backend.np.diag(self.backend.np.diagonal(obs)) + ) + != 0 + ): raise_error(NotImplementedError, "Observable is not diagonal.") keys = list(freq.keys()) if qubit_map is None: @@ -153,7 +158,7 @@ def expectation_from_samples(self, freq, qubit_map=None): for i in qubit_map: index += int(k[qubit_map.index(i)]) * 2 ** (size - 1 - i) expval += obs[index, index] * counts[j] - return np.real(expval) + return self.backend.np.real(expval) def eye(self, dim: Optional[int] = None): """Generate Identity matrix with dimension ``dim``""" @@ -180,8 +185,8 @@ def energy_fluctuation(self, state): energy = self.expectation(state) h = self.matrix h2 = Hamiltonian(nqubits=self.nqubits, matrix=h @ h, backend=self.backend) - average_h2 = h2.expectation(state, normalize=True) - return np.sqrt(np.abs(average_h2 - energy**2)) + average_h2 = self.backend.calculate_expectation_state(h2, state, normalize=True) + return self.backend.np.sqrt(self.backend.np.abs(average_h2 - energy**2)) def __add__(self, o): if isinstance(o, self.__class__): diff --git a/src/qibo/measurements.py b/src/qibo/measurements.py index e46ffe6a00..9f45d27fde 100644 --- a/src/qibo/measurements.py +++ b/src/qibo/measurements.py @@ -168,7 +168,7 @@ def samples(self, binary=True, registers=False): # individual register samples are registered here self.circuit.final_state.samples() if binary: - return self.backend.cast(self._samples, dtype="int32") + return self._samples else: qubits = self.measurement_gate.target_qubits return self.backend.samples_to_decimal(self._samples, len(qubits)) diff --git a/src/qibo/models/dbi/double_bracket.py b/src/qibo/models/dbi/double_bracket.py index e2f8d9b7f2..4c544ea2a7 100644 --- a/src/qibo/models/dbi/double_bracket.py +++ b/src/qibo/models/dbi/double_bracket.py @@ -1,4 +1,4 @@ -from copy import deepcopy +from copy import copy from enum import Enum, auto from typing import Optional @@ -88,7 +88,7 @@ def __init__( ref_state: np.array = None, ): self.h = hamiltonian - self.h0 = deepcopy(self.h) + self.h0 = copy(self.h) self.mode = mode self.scheduling = scheduling self.cost = cost @@ -262,7 +262,7 @@ def loss(self, step: float, d: np.array = None, look_ahead: int = 1): look_ahead (int): number of iteration steps to compute the loss function; """ # copy initial hamiltonian - h_copy = deepcopy(self.h) + h_copy = copy(self.h) for _ in range(look_ahead): self.__call__(mode=self.mode, step=step, d=d) @@ -314,7 +314,9 @@ def cost_expansion(self, d, n): if self.cost is DoubleBracketCostFunction.off_diagonal_norm: coef = off_diagonal_norm_polynomial_expansion_coef(self, d, n) elif self.cost is DoubleBracketCostFunction.least_squares: - coef = least_squares_polynomial_expansion_coef(self, d, n) + coef = least_squares_polynomial_expansion_coef( + self, d, n, backend=self.backend + ) elif self.cost is DoubleBracketCostFunction.energy_fluctuation: coef = energy_fluctuation_polynomial_expansion_coef( self, d, n, self.ref_state diff --git a/src/qibo/models/dbi/utils.py b/src/qibo/models/dbi/utils.py index 10367d71f9..4a04c7a4f8 100644 --- a/src/qibo/models/dbi/utils.py +++ b/src/qibo/models/dbi/utils.py @@ -66,27 +66,28 @@ def str_to_symbolic(name: str): return tensor_op -def cs_angle_sgn(dbi_object, d): +def cs_angle_sgn(dbi_object, d, backend=None): """Calculates the sign of Cauchy-Schwarz Angle :math:`\\langle W(Z), W({\\rm canonical}) \\rangle_{\\rm HS}`.""" - d = dbi_object.backend.cast(d) - norm = np.trace( - np.dot( - np.conjugate( + backend = _check_backend(backend) + d = backend.cast(d) + norm = backend.np.trace( + backend.np.matmul( + backend.np.conj( dbi_object.commutator(dbi_object.diagonal_h_matrix, dbi_object.h.matrix) ).T, dbi_object.commutator(d, dbi_object.h.matrix), ) ) - return np.real(np.sign(norm)) + return backend.np.real(backend.np.sign(norm)) -def decompose_into_pauli_basis(h_matrix: np.array, pauli_operators: list): +def decompose_into_pauli_basis(h_matrix: np.array, pauli_operators: list, backend=None): """finds the decomposition of hamiltonian `h_matrix` into Pauli-Z operators""" nqubits = int(np.log2(h_matrix.shape[0])) - + backend = _check_backend(backend) decomposition = [] for Z_i in pauli_operators: - expect = np.trace(h_matrix @ Z_i) / 2**nqubits + expect = backend.np.trace(h_matrix @ Z_i) / 2**nqubits decomposition.append(expect) return decomposition @@ -184,7 +185,7 @@ def params_to_diagonal_operator( # TODO: write proper tests for normalize=True if normalize: # pragma: no cover d = d / np.linalg.norm(d) - return d + return backend.cast(d) def off_diagonal_norm_polynomial_expansion_coef(dbi_object, d, n): @@ -211,16 +212,18 @@ def off_diagonal_norm_polynomial_expansion_coef(dbi_object, d, n): return coef -def least_squares_polynomial_expansion_coef(dbi_object, d, n: int = 3): +def least_squares_polynomial_expansion_coef(dbi_object, d, n: int = 3, backend=None): """Return the Taylor expansion coefficients of least square cost of `dbi_object.h` and diagonal operator `d` with respect to double bracket rotation duration `s`.""" # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H + backend = _check_backend(backend) Gamma_list = dbi_object.generate_gamma_list(n + 1, d) exp_list = np.array([1 / math.factorial(k) for k in range(n + 1)]) # coefficients coef = np.empty(n) for i in range(n): - coef[i] = np.real( - exp_list[i] * np.trace(dbi_object.backend.cast(d) @ Gamma_list[i + 1]) + coef[i] = backend.np.real( + exp_list[i] + * backend.np.trace(dbi_object.backend.cast(d) @ Gamma_list[i + 1]) ) coef = list(reversed(coef)) return coef @@ -263,3 +266,21 @@ def covariance(a, b): ) coef = list(reversed(coef)) return coef + + +def copy_dbi_object(dbi_object): + """ + Return a copy of the DoubleBracketIteration object. + This is necessary for the `select_best_dbr_generator` function as pytorch do not support deepcopy for leaf tensors. + """ + from copy import copy, deepcopy # pylint: disable=import-outside-toplevel + + dbi_class = dbi_object.__class__ + new = dbi_class.__new__(dbi_class) + + # Manually copy h and h0 as they may be torch tensors + new.h, new.h0 = copy(dbi_object.h), copy(dbi_object.h0) + # Deepcopy the rest of the attributes + for attr in ("mode", "scheduling", "cost", "ref_state"): + setattr(new, attr, deepcopy(getattr(dbi_object, attr, None))) + return new diff --git a/src/qibo/models/dbi/utils_dbr_strategies.py b/src/qibo/models/dbi/utils_dbr_strategies.py index 5aae761fde..fb71cdf49f 100644 --- a/src/qibo/models/dbi/utils_dbr_strategies.py +++ b/src/qibo/models/dbi/utils_dbr_strategies.py @@ -63,9 +63,9 @@ def select_best_dbr_generator( for i, d in enumerate(d_list): # prescribed step durations - dbi_eval = deepcopy(dbi_object) + dbi_eval = copy_dbi_object(dbi_object) d = dbi_eval.backend.cast(d) - flip_list[i] = cs_angle_sgn(dbi_eval, d) + flip_list[i] = cs_angle_sgn(dbi_eval, d, backend=dbi_object.backend) if flip_list[i] != 0: if step is None: step_best = dbi_eval.choose_step( @@ -78,7 +78,7 @@ def select_best_dbr_generator( norms_off_diagonal_restriction[i] = dbi_eval.off_diagonal_norm # canonical if compare_canonical is True: - dbi_eval = deepcopy(dbi_object) + dbi_eval = copy_dbi_object(dbi_object) dbi_eval.mode = DoubleBracketGeneratorType.canonical if step is None: step_best = dbi_eval.choose_step(scheduling=scheduling, **kwargs) @@ -91,7 +91,7 @@ def select_best_dbr_generator( idx_max_loss = np.argmin(norms_off_diagonal_restriction) flip = flip_list[idx_max_loss] step_optimal = optimal_steps[idx_max_loss] - dbi_eval = deepcopy(dbi_object) + dbi_eval = copy_dbi_object(dbi_object) if idx_max_loss == len(d_list) and compare_canonical is True: # canonical dbi_eval(step=step_optimal, mode=DoubleBracketGeneratorType.canonical) @@ -126,13 +126,17 @@ def gradient_numerical( nqubits = dbi_object.nqubits grad = np.zeros(len(d_params)) d = params_to_diagonal_operator( - d_params, nqubits, parameterization=parameterization, **kwargs + d_params, nqubits, parameterization=parameterization, **kwargs, backend=backend ) for i in range(len(d_params)): - params_new = backend.cast(d_params, copy=True) - params_new[i] += delta + params_new = backend.to_numpy(d_params).copy() + params_new[i] = params_new[i] + delta d_new = params_to_diagonal_operator( - params_new, nqubits, parameterization=parameterization, **kwargs + params_new, + nqubits, + parameterization=parameterization, + **kwargs, + backend=backend, ) # find the increment of a very small step grad[i] = (dbi_object.loss(s, d_new) - dbi_object.loss(s, d)) / delta @@ -230,6 +234,7 @@ def gradient_descent( parameterization=parameterization, pauli_operator_dict=pauli_operator_dict, normalize=normalize, + backend=backend, ) loss_hist = [dbi_object.loss(0.0, d=d)] d_params_hist = [d_params] @@ -266,6 +271,7 @@ def func_loss_to_lr(lr): parameterization=parameterization, pauli_operator_dict=pauli_operator_dict, normalize=normalize, + backend=backend, ) return dbi_object.loss(step=s, d=d_eval) @@ -288,6 +294,7 @@ def func_loss_to_lr(lr): parameterization=parameterization, pauli_operator_dict=pauli_operator_dict, normalize=normalize, + backend=backend, ) s = dbi_object.choose_step(d=d) dbi_object(step=s, d=d) diff --git a/src/qibo/models/encodings.py b/src/qibo/models/encodings.py index 611ffd9534..3dd0b88dd7 100644 --- a/src/qibo/models/encodings.py +++ b/src/qibo/models/encodings.py @@ -8,6 +8,7 @@ from qibo import gates from qibo.config import raise_error +from qibo.gates.gates import _check_engine from qibo.models.circuit import Circuit @@ -435,8 +436,9 @@ def _generate_rbs_angles(data, nqubits: int, architecture: str): list: list of phases for RBS gates. """ if architecture == "diagonal": + engine = _check_engine(data) phases = [ - math.atan2(np.linalg.norm(data[k + 1 :]), data[k]) + math.atan2(engine.linalg.norm(data[k + 1 :]), data[k]) for k in range(len(data) - 2) ] phases.append(math.atan2(data[-1], data[-2])) diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py index 62e52bb0a2..314b946938 100644 --- a/src/qibo/models/error_mitigation.py +++ b/src/qibo/models/error_mitigation.py @@ -229,17 +229,20 @@ def sample_training_circuit_cdr( replacement.append(rep_gates) distance.append( - np.linalg.norm( - gate.matrix(backend) - - backend.cast([rep_gate.matrix(backend) for rep_gate in rep_gates]), - ord="fro", - axis=(1, 2), + backend.np.real( + backend.np.linalg.norm( + gate.matrix(backend) + - backend.cast( + [rep_gate.matrix(backend) for rep_gate in rep_gates] + ), + ord="fro", + axis=(1, 2), + ) ) ) - distance = np.vstack(distance) - prob = np.exp(-(distance**2) / sigma**2) - prob = backend.cast(prob, dtype=prob.dtype) + distance = backend.np.vstack(distance) + prob = backend.np.exp(-(distance**2) / sigma**2) index = local_state.choice( range(len(gates_to_replace)), diff --git a/src/qibo/models/variational.py b/src/qibo/models/variational.py index aeeaa6352c..7bdc828bac 100644 --- a/src/qibo/models/variational.py +++ b/src/qibo/models/variational.py @@ -93,7 +93,11 @@ def minimize( # TODO: check if we can use this shortcut # dtype = getattr(self.hamiltonian.backend.np, self.hamiltonian.backend._dtypes.get('DTYPE')) dtype = self.hamiltonian.backend.np.float64 - loss = lambda p, c, h: dtype(loss_func(p, c, h)) + loss = ( + (lambda p, c, h: loss_func(p, c, h).item()) + if str(dtype) == "torch.float64" + else (lambda p, c, h: dtype(loss_func(p, c, h))) + ) elif method != "sgd": loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(loss_func(p, c, h)) result, parameters, extra = self.optimizers.optimize( diff --git a/src/qibo/optimizers.py b/src/qibo/optimizers.py index 6645326d95..2696deda2b 100644 --- a/src/qibo/optimizers.py +++ b/src/qibo/optimizers.py @@ -241,7 +241,9 @@ def sgd( """Stochastic Gradient Descent (SGD) optimizer using Tensorflow backpropagation. See `tf.keras.Optimizers `_ - for a list of the available optimizers. + for a list of the available optimizers for Tensorflow. + See `torch.optim `_ for a list of the available + optimizers for PyTorch. Args: loss (callable): Loss as a function of variational parameters to be @@ -260,9 +262,6 @@ def sgd( a message of the loss function. """ - if not backend.name == "tensorflow": - raise_error(RuntimeError, "SGD optimizer requires Tensorflow backend.") - sgd_options = { "nepochs": 1000000, "nmessage": 1000, @@ -272,7 +271,53 @@ def sgd( if options is not None: sgd_options.update(options) - # proceed with the training + if backend.name == "tensorflow": + return _sgd_tf( + loss, + initial_parameters, + args, + sgd_options, + compile, + backend, + callback=callback, + ) + + if backend.name == "pytorch": + if compile: + log.warning( + "PyTorch does not support compilation of the optimization graph." + ) + return _sgd_torch( + loss, initial_parameters, args, sgd_options, backend, callback=callback + ) + + raise_error(RuntimeError, "SGD optimizer requires Tensorflow or PyTorch backend.") + + +def _sgd_torch(loss, initial_parameters, args, sgd_options, backend, callback=None): + + vparams = initial_parameters + optimizer = getattr(backend.np.optim, sgd_options["optimizer"])( + params=[vparams], lr=sgd_options["learning_rate"] + ) + + for e in range(sgd_options["nepochs"]): + optimizer.zero_grad() + l = loss(vparams, *args) + l.backward() + optimizer.step() + if callback is not None: + callback(backend.to_numpy(vparams)) + if e % sgd_options["nmessage"] == 1: + log.info("ite %d : loss %f", e, l.item()) + + return loss(vparams, *args).item(), vparams.detach().numpy(), sgd_options + + +def _sgd_tf( + loss, initial_parameters, args, sgd_options, compile, backend, callback=None +): + vparams = backend.tf.Variable(initial_parameters) optimizer = getattr(backend.tf.optimizers, sgd_options["optimizer"])( learning_rate=sgd_options["learning_rate"] diff --git a/src/qibo/quantum_info/basis.py b/src/qibo/quantum_info/basis.py index 24b2a17507..5ea541bca0 100644 --- a/src/qibo/quantum_info/basis.py +++ b/src/qibo/quantum_info/basis.py @@ -100,28 +100,33 @@ def pauli_basis( else: basis_full = basis_single + basis_full = backend.cast(basis_full, dtype=basis_full[0].dtype) + if vectorize and sparse: basis, indexes = [], [] for row in basis_full: row = vectorization(row, order=order, backend=backend) - row_indexes = list(np.flatnonzero(row)) + row_indexes = backend.np.flatnonzero(row) indexes.append(row_indexes) basis.append(row[row_indexes]) del row elif vectorize and not sparse: basis = [ - vectorization(matrix, order=order, backend=backend) for matrix in basis_full + vectorization( + backend.cast(matrix, dtype=matrix.dtype), order=order, backend=backend + ) + for matrix in basis_full ] else: basis = basis_full - basis = backend.cast(basis, dtype=backend.dtype) + basis = backend.cast(basis, dtype=basis[0].dtype) if normalize: - basis /= np.sqrt(2**nqubits) + basis = basis / np.sqrt(2**nqubits) if vectorize and sparse: - indexes = backend.cast(indexes) + indexes = backend.cast(indexes, dtype=indexes[0][0].dtype) return basis, indexes @@ -198,7 +203,7 @@ def comp_basis_to_pauli( pauli_order=pauli_order, backend=backend, ) - elements = np.conj(elements) + elements = backend.np.conj(elements) return elements, indexes @@ -212,8 +217,7 @@ def comp_basis_to_pauli( backend=backend, ) - unitary = np.conj(unitary) - unitary = backend.cast(unitary, dtype=unitary.dtype) + unitary = backend.np.conj(unitary) return unitary @@ -267,12 +271,12 @@ def pauli_to_comp_basis( pauli_order=pauli_order, backend=backend, ) - unitary = np.transpose(unitary) + unitary = unitary.T if sparse: elements, indexes = [], [] for row in unitary: - index_list = list(np.flatnonzero(row)) + index_list = backend.np.flatnonzero(row) indexes.append(index_list) elements.append(row[index_list]) diff --git a/src/qibo/quantum_info/entanglement.py b/src/qibo/quantum_info/entanglement.py index 28cbb70305..e1f124bd64 100644 --- a/src/qibo/quantum_info/entanglement.py +++ b/src/qibo/quantum_info/entanglement.py @@ -50,7 +50,7 @@ def concurrence(state, bipartition, check_purity: bool = True, backend=None): nqubits = int(np.log2(state.shape[0])) if check_purity is True: - purity_total_system = purity(state) + purity_total_system = purity(state, backend=backend) mixed = bool(abs(purity_total_system - 1.0) > PRECISION_TOL) if mixed is True: @@ -65,7 +65,7 @@ def concurrence(state, bipartition, check_purity: bool = True, backend=None): else backend.partial_trace_density_matrix(state, bipartition, nqubits) ) - purity_reduced = purity(reduced_density_matrix) + purity_reduced = purity(reduced_density_matrix, backend=backend) if purity_reduced - 1.0 > 0.0: purity_reduced = round(purity_reduced, 7) @@ -229,7 +229,7 @@ def meyer_wallach_entanglement(circuit, backend=None): rho_r = backend.partial_trace_density_matrix(rho, trace_q, nqubits) - trace = purity(rho_r) + trace = purity(rho_r, backend=backend) ent += trace diff --git a/src/qibo/quantum_info/entropies.py b/src/qibo/quantum_info/entropies.py index d6499e15dc..20ff3e2af2 100644 --- a/src/qibo/quantum_info/entropies.py +++ b/src/qibo/quantum_info/entropies.py @@ -6,6 +6,7 @@ from scipy.linalg import fractional_matrix_power from qibo.backends import _check_backend +from qibo.backends.pytorch import PyTorchBackend from qibo.config import PRECISION_TOL, raise_error from qibo.quantum_info.metrics import _check_hermitian_or_not_gpu, purity @@ -56,17 +57,19 @@ def shannon_entropy(prob_dist, base: float = 2, backend=None): total_sum = backend.np.sum(prob_dist) - if np.abs(total_sum - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum) - 1.0) > PRECISION_TOL: raise_error(ValueError, "Probability array must sum to 1.") - log_prob = np.where(prob_dist != 0, np.log2(prob_dist) / np.log2(base), 0.0) + log_prob = backend.np.where( + prob_dist != 0, backend.np.log2(prob_dist) / np.log2(base), 0.0 + ) shan_entropy = -backend.np.sum(prob_dist * log_prob) # absolute value if entropy == 0.0 to avoid returning -0.0 - shan_entropy = np.abs(shan_entropy) if shan_entropy == 0.0 else shan_entropy + shan_entropy = backend.np.abs(shan_entropy) if shan_entropy == 0.0 else shan_entropy - return complex(shan_entropy).real + return np.real(float(shan_entropy)) def classical_relative_entropy(prob_dist_p, prob_dist_q, base: float = 2, backend=None): @@ -93,13 +96,8 @@ def classical_relative_entropy(prob_dist_p, prob_dist_q, base: float = 2, backen float: Classical relative entropy between :math:`\\mathbf{p}` and :math:`\\mathbf{q}`. """ backend = _check_backend(backend) - - if isinstance(prob_dist_p, list): - # np.float64 is necessary instead of native float because of tensorflow - prob_dist_p = backend.cast(prob_dist_p, dtype=np.float64) - if isinstance(prob_dist_q, list): - # np.float64 is necessary instead of native float because of tensorflow - prob_dist_q = backend.cast(prob_dist_q, dtype=np.float64) + prob_dist_p = backend.cast(prob_dist_p, dtype=np.float64) + prob_dist_q = backend.cast(prob_dist_q, dtype=np.float64) if (len(prob_dist_p.shape) != 1) or (len(prob_dist_q.shape) != 1): raise_error( @@ -125,19 +123,19 @@ def classical_relative_entropy(prob_dist_p, prob_dist_q, base: float = 2, backen total_sum_q = backend.np.sum(prob_dist_q) - if np.abs(total_sum_p - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum_p) - 1.0) > PRECISION_TOL: raise_error(ValueError, "First probability array must sum to 1.") - if np.abs(total_sum_q - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum_q) - 1.0) > PRECISION_TOL: raise_error(ValueError, "Second probability array must sum to 1.") entropy_p = -1 * shannon_entropy(prob_dist_p, base=base, backend=backend) - log_prob_q = np.where( - prob_dist_q != 0.0, np.log2(prob_dist_q) / np.log2(base), -np.inf + log_prob_q = backend.np.where( + prob_dist_q != 0.0, backend.np.log2(prob_dist_q) / np.log2(base), -np.inf ) - log_prob = np.where(prob_dist_p != 0.0, log_prob_q, 0.0) + log_prob = backend.np.where(prob_dist_p != 0.0, log_prob_q, 0.0) relative = backend.np.sum(prob_dist_p * log_prob) @@ -181,10 +179,7 @@ def classical_renyi_entropy( float: Classical Rényi entropy :math:`H_{\\alpha}`. """ backend = _check_backend(backend) - - if isinstance(prob_dist, list): - # np.float64 is necessary instead of native float because of tensorflow - prob_dist = backend.cast(prob_dist, dtype=np.float64) + prob_dist = backend.cast(prob_dist, dtype=np.float64) if not isinstance(alpha, (float, int)): raise_error( @@ -214,7 +209,7 @@ def classical_renyi_entropy( total_sum = backend.np.sum(prob_dist) - if np.abs(total_sum - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum) - 1.0) > PRECISION_TOL: raise_error(ValueError, "Probability array must sum to 1.") if alpha == 0.0: @@ -224,11 +219,11 @@ def classical_renyi_entropy( return shannon_entropy(prob_dist, base=base, backend=backend) if alpha == np.inf: - return -1 * np.log2(max(prob_dist)) / np.log2(base) + return -1 * backend.np.log2(max(prob_dist)) / np.log2(base) total_sum = backend.np.sum(prob_dist**alpha) - renyi_ent = (1 / (1 - alpha)) * np.log2(total_sum) / np.log2(base) + renyi_ent = (1 / (1 - alpha)) * backend.np.log2(total_sum) / np.log2(base) return renyi_ent @@ -273,13 +268,8 @@ def classical_relative_renyi_entropy( float: Classical relative Rényi entropy :math:`H_{\\alpha}(\\mathbf{p} \\, \\| \\, \\mathbf{q})`. """ backend = _check_backend(backend) - - if isinstance(prob_dist_p, list): - # np.float64 is necessary instead of native float because of tensorflow - prob_dist_p = backend.cast(prob_dist_p, dtype=np.float64) - if isinstance(prob_dist_q, list): - # np.float64 is necessary instead of native float because of tensorflow - prob_dist_q = backend.cast(prob_dist_q, dtype=np.float64) + prob_dist_p = backend.cast(prob_dist_p, dtype=np.float64) + prob_dist_q = backend.cast(prob_dist_q, dtype=np.float64) if (len(prob_dist_p.shape) != 1) or (len(prob_dist_q.shape) != 1): raise_error( @@ -313,17 +303,17 @@ def classical_relative_renyi_entropy( total_sum_p = backend.np.sum(prob_dist_p) total_sum_q = backend.np.sum(prob_dist_q) - if np.abs(total_sum_p - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum_p) - 1.0) > PRECISION_TOL: raise_error(ValueError, "First probability array must sum to 1.") - if np.abs(total_sum_q - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum_q) - 1.0) > PRECISION_TOL: raise_error(ValueError, "Second probability array must sum to 1.") if alpha == 0.5: - total_sum = np.sqrt(prob_dist_p * prob_dist_q) + total_sum = backend.np.sqrt(prob_dist_p * prob_dist_q) total_sum = backend.np.sum(total_sum) - return -2 * np.log2(total_sum) / np.log2(base) + return -2 * backend.np.log2(total_sum) / np.log2(base) if alpha == 1.0: return classical_relative_entropy( @@ -331,14 +321,14 @@ def classical_relative_renyi_entropy( ) if alpha == np.inf: - return np.log2(max(prob_dist_p / prob_dist_q)) / np.log2(base) + return backend.np.log2(max(prob_dist_p / prob_dist_q)) / np.log2(base) prob_p = prob_dist_p**alpha prob_q = prob_dist_q ** (1 - alpha) total_sum = backend.np.sum(prob_p * prob_q) - return (1 / (alpha - 1)) * np.log2(total_sum) / np.log2(base) + return (1 / (alpha - 1)) * backend.np.log2(total_sum) / np.log2(base) def classical_tsallis_entropy(prob_dist, alpha: float, base: float = 2, backend=None): @@ -396,7 +386,7 @@ def classical_tsallis_entropy(prob_dist, alpha: float, base: float = 2, backend= total_sum = backend.np.sum(prob_dist) - if np.abs(total_sum - 1.0) > PRECISION_TOL: + if np.abs(float(total_sum) - 1.0) > PRECISION_TOL: raise_error(ValueError, "Probability array must sum to 1.") if alpha == 1.0: @@ -459,24 +449,28 @@ def von_neumann_entropy( f"check_hermitian must be type bool, but it is type {type(check_hermitian)}.", ) - if purity(state) == 1.0: + if purity(state, backend=backend) == 1.0: if return_spectrum: return 0.0, backend.cast([0.0], dtype=float) return 0.0 - if not check_hermitian or _check_hermitian_or_not_gpu(state, backend=backend): - eigenvalues = np.linalg.eigvalsh(state) - else: - eigenvalues = np.linalg.eigvals(state) + eigenvalues = backend.calculate_eigenvalues( + state, + hermitian=( + not check_hermitian or _check_hermitian_or_not_gpu(state, backend=backend) + ), + ) - log_prob = np.where(eigenvalues > 0, np.log2(eigenvalues) / np.log2(base), 0.0) + log_prob = backend.np.where( + backend.np.real(eigenvalues) > 0.0, + backend.np.log2(eigenvalues) / np.log2(base), + 0.0, + ) - ent = -np.sum(eigenvalues * log_prob) + ent = -backend.np.sum(eigenvalues * log_prob) # absolute value if entropy == 0.0 to avoid returning -0.0 - ent = np.abs(ent) if ent == 0.0 else ent - - ent = float(ent) + ent = backend.np.abs(ent) if ent == 0.0 else backend.np.real(ent) if return_spectrum: log_prob = backend.cast(log_prob, dtype=log_prob.dtype) @@ -511,6 +505,8 @@ def relative_von_neumann_entropy( float: Relative (von-Neumann) entropy :math:`S(\\rho \\, \\| \\, \\sigma)`. """ backend = _check_backend(backend) + state = backend.cast(state) + target = backend.cast(target) if ( (len(state.shape) >= 3) @@ -541,39 +537,46 @@ def relative_von_neumann_entropy( f"check_hermitian must be type bool, but it is type {type(check_hermitian)}.", ) - if purity(state) == 1.0 and purity(target) == 1.0: + if purity(state, backend=backend) == 1.0 and purity(target, backend=backend) == 1.0: return 0.0 if len(state.shape) == 1: - state = np.outer(state, np.conj(state)) + state = backend.np.outer(state, backend.np.conj(state)) if len(target.shape) == 1: - target = np.outer(target, np.conj(target)) + target = backend.np.outer(target, backend.np.conj(target)) - if not check_hermitian or _check_hermitian_or_not_gpu(state, backend=backend): - eigenvalues_state = np.linalg.eigvalsh(state) - else: - eigenvalues_state = np.linalg.eigvals(state) - - if not check_hermitian or _check_hermitian_or_not_gpu(target, backend=backend): - eigenvalues_target = np.linalg.eigvalsh(target) - else: - eigenvalues_target = np.linalg.eigvals(target) + eigenvalues_state = backend.calculate_eigenvalues( + state, + hermitian=( + not check_hermitian or _check_hermitian_or_not_gpu(state, backend=backend) + ), + ) + eigenvalues_target = backend.calculate_eigenvalues( + target, + hermitian=( + not check_hermitian or _check_hermitian_or_not_gpu(target, backend=backend) + ), + ) - log_state = np.where( - eigenvalues_state > 0, np.log2(eigenvalues_state) / np.log2(base), 0.0 + log_state = backend.np.where( + backend.np.real(eigenvalues_state) > 0, + backend.np.log2(eigenvalues_state) / np.log2(base), + 0.0, ) - log_target = np.where( - eigenvalues_target > 0, np.log2(eigenvalues_target) / np.log2(base), -np.inf + log_target = backend.np.where( + backend.np.real(eigenvalues_target) > 0, + backend.np.log2(eigenvalues_target) / np.log2(base), + -np.inf, ) - log_target = np.where(eigenvalues_state != 0.0, log_target, 0.0) + log_target = backend.np.where(eigenvalues_state != 0.0, log_target, 0.0) - entropy_state = np.sum(eigenvalues_state * log_state) + entropy_state = backend.np.sum(eigenvalues_state * log_state) - relative = np.sum(eigenvalues_state * log_target) + relative = backend.np.sum(eigenvalues_state * log_target) - return float(entropy_state - relative) + return float(backend.np.real(entropy_state - relative)) def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None): @@ -632,7 +635,7 @@ def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None if base <= 0.0: raise_error(ValueError, "log base must be non-negative.") - if abs(purity(state) - 1.0) < PRECISION_TOL: + if abs(purity(state, backend=backend) - 1.0) < PRECISION_TOL: return 0.0 if alpha == 0.0: @@ -644,11 +647,11 @@ def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None if alpha == np.inf: return ( -1 - * np.log2(backend.calculate_norm_density_matrix(state, order=2)) + * backend.np.log2(backend.calculate_norm_density_matrix(state, order=2)) / np.log2(base) ) - log = np.log2(np.trace(_matrix_power(state, alpha, backend))) + log = backend.np.log2(backend.np.trace(_matrix_power(state, alpha, backend))) return (1 / (1 - alpha)) * log / np.log2(base) @@ -694,7 +697,8 @@ def relative_renyi_entropy( float: Relative Rényi entropy :math:`H_{\\alpha}(\\rho \\, \\| \\, \\sigma)`. """ backend = _check_backend(backend) - + state = backend.cast(state) + target = backend.cast(target) if ( (len(state.shape) >= 3) or (len(state) == 0) @@ -726,9 +730,9 @@ def relative_renyi_entropy( if base <= 0.0: raise_error(ValueError, "log base must be non-negative.") - purity_target = purity(target) + purity_target = purity(target, backend=backend) if ( - abs(purity(state) - 1.0) < PRECISION_TOL + abs(purity(state, backend=backend) - 1.0) < PRECISION_TOL and abs(purity_target - 1) < PRECISION_TOL ): return 0.0 @@ -740,7 +744,7 @@ def relative_renyi_entropy( ) if len(state.shape) == 1: - state = np.outer(state, np.conj(state)) + state = backend.np.outer(state, backend.np.conj(state)) if alpha == 1.0: return relative_von_neumann_entropy(state, target, base, backend=backend) @@ -749,7 +753,7 @@ def relative_renyi_entropy( new_state = _matrix_power(state, 0.5, backend) new_target = _matrix_power(target, 0.5, backend) - log = np.log2( + log = backend.np.log2( backend.calculate_norm_density_matrix(new_state @ new_target, order=1) ) @@ -757,7 +761,7 @@ def relative_renyi_entropy( log = _matrix_power(state, alpha, backend) log = log @ _matrix_power(target, 1 - alpha, backend) - log = np.log2(np.trace(log)) + log = backend.np.log2(backend.np.trace(log)) return (1 / (alpha - 1)) * log / np.log2(base) @@ -807,13 +811,15 @@ def tsallis_entropy(state, alpha: float, base: float = 2, backend=None): if base <= 0.0: raise_error(ValueError, "log base must be non-negative.") - if abs(purity(state) - 1.0) < PRECISION_TOL: + if abs(purity(state, backend=backend) - 1.0) < PRECISION_TOL: return 0.0 if alpha == 1.0: return von_neumann_entropy(state, base=base, backend=backend) - return (1 / (1 - alpha)) * (np.trace(_matrix_power(state, alpha, backend)) - 1) + return (1 / (1 - alpha)) * ( + backend.np.trace(_matrix_power(state, alpha, backend)) - 1 + ) def entanglement_entropy( @@ -890,11 +896,9 @@ def _matrix_power(matrix, alpha, backend): ]: # pragma: no cover new_matrix = backend.to_numpy(matrix) else: - new_matrix = np.copy(matrix) + new_matrix = backend.np.copy(matrix) if len(new_matrix.shape) == 1: - new_matrix = np.outer(new_matrix, np.conj(new_matrix)) - - new_matrix = fractional_matrix_power(new_matrix, alpha) + new_matrix = backend.np.outer(new_matrix, backend.np.conj(new_matrix)) - return backend.cast(new_matrix, dtype=new_matrix.dtype) + return backend.cast(fractional_matrix_power(backend.to_numpy(new_matrix), alpha)) diff --git a/src/qibo/quantum_info/metrics.py b/src/qibo/quantum_info/metrics.py index c7f5f91a16..4699efbda6 100644 --- a/src/qibo/quantum_info/metrics.py +++ b/src/qibo/quantum_info/metrics.py @@ -9,7 +9,7 @@ from qibo.config import PRECISION_TOL, raise_error -def purity(state): +def purity(state, backend=None): """Purity of a quantum state :math:`\\rho`. This is given by @@ -22,7 +22,7 @@ def purity(state): Returns: float: Purity of quantum ``state`` :math:`\\rho`. """ - + backend = _check_backend(backend) if ( (len(state.shape) >= 3) or (len(state) == 0) @@ -34,18 +34,13 @@ def purity(state): ) if len(state.shape) == 1: - pur = np.abs(np.dot(np.conj(state), state)) ** 2 + pur = backend.np.real(backend.calculate_norm(state)) ** 2 else: - pur = np.real(np.trace(np.dot(state, state))) - - # this is necessary to remove the float from inside - # a 0-dim ndarray - pur = float(pur) + pur = backend.np.real(backend.np.trace(backend.np.matmul(state, state))) + return float(pur) - return pur - -def impurity(state): +def impurity(state, backend=None): """Impurity of quantum state :math:`\\rho`. This is given by :math:`1 - \\text{purity}(\\rho)`, where :math:`\\text{purity}` @@ -57,7 +52,7 @@ def impurity(state): Returns: float: impurity of ``state`` :math:`\\rho`. """ - return 1 - purity(state) + return 1 - purity(state, backend=backend) def trace_distance(state, target, check_hermitian: bool = False, backend=None): @@ -107,15 +102,17 @@ def trace_distance(state, target, check_hermitian: bool = False, backend=None): ) if len(state.shape) == 1: - state = np.outer(np.conj(state), state) - target = np.outer(np.conj(target), target) + state = backend.np.outer(backend.np.conj(state), state) + target = backend.np.outer(backend.np.conj(target), target) difference = state - target if check_hermitian is True: hermitian = bool( float( backend.calculate_norm_density_matrix( - np.transpose(np.conj(difference)) - difference, order=2 + backend.np.transpose(backend.np.conj(difference), (1, 0)) + - difference, + order=2, ) ) <= PRECISION_TOL @@ -128,21 +125,16 @@ def trace_distance(state, target, check_hermitian: bool = False, backend=None): "CupyBackend does not support `np.linalg.eigvals`" + "for non-Hermitian `state - target`.", ) - eigenvalues = ( - np.linalg.eigvalsh(difference) - if hermitian - else np.linalg.eigvals(difference) - ) + eigenvalues = backend.calculate_eigenvalues(difference, hermitian=hermitian) else: - eigenvalues = np.linalg.eigvalsh(difference) + eigenvalues = backend.calculate_eigenvalues(difference) - distance = np.sum(np.absolute(eigenvalues)) / 2 - distance = float(distance) + distance = backend.np.sum(backend.np.absolute(eigenvalues)) / 2 return distance -def hilbert_schmidt_distance(state, target): +def hilbert_schmidt_distance(state, target, backend=None): """Hilbert-Schmidt distance between two quantum states: .. math:: @@ -152,11 +144,15 @@ def hilbert_schmidt_distance(state, target): Args: state (ndarray): statevector or density matrix. target (ndarray): statevector or density matrix. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Returns: float: Hilbert-Schmidt distance between ``state`` :math:`\\rho` and ``target`` :math:`\\sigma`. """ + backend = _check_backend(backend) if state.shape != target.shape: raise_error( @@ -172,10 +168,10 @@ def hilbert_schmidt_distance(state, target): ) if len(state.shape) == 1: - state = np.outer(np.conj(state), state) - target = np.outer(np.conj(target), target) + state = backend.np.outer(backend.np.conj(state), state) + target = backend.np.outer(backend.np.conj(target), target) - distance = np.real(np.trace((state - target) ** 2)) + distance = backend.np.real(backend.np.trace((state - target) ** 2)) distance = float(distance) return distance @@ -207,7 +203,8 @@ def fidelity(state, target, check_hermitian: bool = False, backend=None): float: Fidelity between ``state`` :math:`\\rho` and ``target`` :math:`\\sigma`. """ backend = _check_backend(backend) - + state = backend.cast(state, dtype=state.dtype) + target = backend.cast(target, dtype=target.dtype) if state.shape != target.shape: raise_error( TypeError, @@ -229,60 +226,63 @@ def fidelity(state, target, check_hermitian: bool = False, backend=None): # check purity if both states are density matrices if len(state.shape) == 2 and len(target.shape) == 2: - purity_state = purity(state) - purity_target = purity(target) + purity_state = purity(state, backend=backend) + purity_target = purity(target, backend=backend) # if both states are mixed, default to full fidelity calculation if ( abs(purity_state - 1) > PRECISION_TOL and abs(purity_target - 1) > PRECISION_TOL ): - # using eigh since rho is supposed to be Hermitian - if check_hermitian is False or _check_hermitian_or_not_gpu( + hermitian = check_hermitian is False or _check_hermitian_or_not_gpu( state, backend=backend - ): - eigenvalues, eigenvectors = np.linalg.eigh(state) - else: - eigenvalues, eigenvectors = np.linalg.eig(state) + ) + # using eigh since rho is supposed to be Hermitian + eigenvalues, eigenvectors = backend.calculate_eigenvectors( + state, hermitian=hermitian + ) state = np.zeros(state.shape, dtype=complex) state = backend.cast(state, dtype=state.dtype) - for eig, eigvec in zip(eigenvalues, np.transpose(eigenvectors)): - matrix = np.sqrt(eig) * np.outer(eigvec, np.conj(eigvec)) + for eig, eigvec in zip( + eigenvalues, backend.np.transpose(eigenvectors, (1, 0)) + ): + matrix = backend.np.sqrt(eig) * backend.np.outer( + eigvec, backend.np.conj(eigvec) + ) matrix = backend.cast(matrix, dtype=matrix.dtype) - state += matrix + state = state + matrix del matrix fid = state @ target @ state # since sqrt(rho) is Hermitian, we can use eigh again - if check_hermitian is False or _check_hermitian_or_not_gpu( - fid, backend=backend - ): - eigenvalues, eigenvectors = np.linalg.eigh(fid) - else: - eigenvalues, eigenvectors = np.linalg.eig(fid) + eigenvalues, eigenvectors = backend.calculate_eigenvectors( + fid, hermitian=hermitian + ) fid = np.zeros(state.shape, dtype=complex) fid = backend.cast(fid, dtype=fid.dtype) - for eig, eigvec in zip(eigenvalues, np.transpose(eigenvectors)): - if eig > PRECISION_TOL: - matrix = np.sqrt(eig) * np.outer(eigvec, np.conj(eigvec)) + for eig, eigvec in zip( + eigenvalues, backend.np.transpose(eigenvectors, (1, 0)) + ): + if backend.np.real(eig) > PRECISION_TOL: + matrix = backend.np.sqrt(eig) * backend.np.outer( + eigvec, backend.np.conj(eigvec) + ) matrix = backend.cast(matrix, dtype=matrix.dtype) - fid += matrix + fid = fid + matrix del matrix - fid = np.real(np.trace(fid)) ** 2 + fid = backend.np.real(backend.np.trace(fid)) ** 2 return fid # if any of the states is pure, perform lighter calculation fid = ( - np.abs(np.dot(np.conj(state), target)) ** 2 + backend.np.abs(backend.np.matmul(backend.np.conj(state), target)) ** 2 if len(state.shape) == 1 - else np.real(np.trace(np.dot(state, target))) + else backend.np.real(backend.np.trace(backend.np.matmul(state, target))) ) - fid = float(fid) - return fid @@ -333,8 +333,10 @@ def bures_angle(state, target, check_hermitian: bool = False, backend=None): Returns: float: Bures angle between ``state`` and ``target``. """ - angle = np.arccos( - np.sqrt(fidelity(state, target, check_hermitian, backend=backend)) + backend = _check_backend(backend) + + angle = backend.np.arccos( + backend.np.sqrt(fidelity(state, target, check_hermitian, backend=backend)) ) return angle @@ -362,9 +364,11 @@ def bures_distance(state, target, check_hermitian: bool = False, backend=None): Returns: float: Bures distance between ``state`` and ``target``. """ - distance = np.sqrt( - 2 * (1 - np.sqrt(fidelity(state, target, check_hermitian, backend=backend))) + backend = _check_backend(backend) + sqrt_fid = backend.np.sqrt( + fidelity(state, target, check_hermitian, backend=backend) ) + distance = backend.np.sqrt(2 * (1 - sqrt_fid)) return distance @@ -404,7 +408,10 @@ def process_fidelity(channel, target=None, check_unitary: bool = False, backend= if check_unitary is True: norm_channel = float( backend.calculate_norm_density_matrix( - np.dot(np.conj(np.transpose(channel)), channel) - np.eye(dim**2) + backend.np.matmul( + backend.np.conj(backend.np.transpose(channel, (1, 0))), channel + ) + - backend.np.eye(dim**2) ) ) if target is None and norm_channel > PRECISION_TOL: @@ -412,7 +419,10 @@ def process_fidelity(channel, target=None, check_unitary: bool = False, backend= if target is not None: norm_target = float( backend.calculate_norm( - np.dot(np.conj(np.transpose(target)), target) - np.eye(dim**2) + backend.np.matmul( + backend.np.conj(backend.np.transpose(target, (1, 0))), target + ) + - backend.np.eye(dim**2) ) ) if (norm_channel > PRECISION_TOL) and (norm_target > PRECISION_TOL): @@ -420,14 +430,15 @@ def process_fidelity(channel, target=None, check_unitary: bool = False, backend= if target is None: # With no target, return process fidelity with Identity channel - process_fid = np.real(np.trace(channel)) / dim**2 + process_fid = backend.np.real(backend.np.trace(channel)) / dim**2 process_fid = float(process_fid) return process_fid - process_fid = np.dot(np.transpose(np.conj(channel)), target) - process_fid = np.real(np.trace(process_fid)) / dim**2 - process_fid = float(process_fid) + process_fid = backend.np.matmul( + backend.np.transpose(backend.np.conj(channel), (1, 0)), target + ) + process_fid = backend.np.real(backend.np.trace(process_fid)) / dim**2 return process_fid @@ -567,6 +578,9 @@ def diamond_norm(channel, target=None, backend=None, **kwargs): # pragma: no co channel (ndarray): row-vectorized Choi representation of a quantum channel. target (ndarray, optional): row-vectorized Choi representation of a target quantum channel. Defaults to ``None``. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. kwargs: optional arguments to pass to CVXPY solver. For more information, please visit `CVXPY's API documentation `_. @@ -579,6 +593,8 @@ def diamond_norm(channel, target=None, backend=None, **kwargs): # pragma: no co """ import cvxpy + backend = _check_backend(backend) + if target is not None: if channel.shape != target.shape: raise_error( @@ -595,9 +611,9 @@ def diamond_norm(channel, target=None, backend=None, **kwargs): # pragma: no co channel = backend.to_numpy(channel) - channel = np.transpose(channel) - channel_real = np.real(channel) - channel_imag = np.imag(channel) + channel = backend.np.transpose(channel, (1, 0)) + channel_real = backend.np.real(channel) + channel_imag = backend.np.imag(channel) dim = int(np.sqrt(channel.shape[0])) @@ -787,8 +803,10 @@ def frame_potential( unitary_2.set_parameters(params_2) unitary_2 = unitary_2.unitary(backend) / np.sqrt(dim) - potential += np.abs( - np.trace(np.transpose(np.conj(unitary_1)) @ unitary_2) + potential += backend.np.abs( + backend.np.trace( + backend.np.transpose(backend.np.conj(unitary_1), (1, 0)) @ unitary_2 + ) ) ** (2 * power_t) return potential / samples**2 @@ -815,7 +833,7 @@ def _check_hermitian_or_not_gpu(matrix, backend=None): backend = _check_backend(backend) norm = backend.calculate_norm_density_matrix( - np.transpose(np.conj(matrix)) - matrix, order=2 + backend.np.conj(matrix).T - matrix, order=2 ) hermitian = bool(float(norm) <= PRECISION_TOL) diff --git a/src/qibo/quantum_info/quantum_networks.py b/src/qibo/quantum_info/quantum_networks.py index e7bf52e0ad..48b34a0ff2 100644 --- a/src/qibo/quantum_info/quantum_networks.py +++ b/src/qibo/quantum_info/quantum_networks.py @@ -119,10 +119,11 @@ def is_hermitian( self._pure = False reshaped = self._backend.cast( - np.reshape(self._matrix, (self.dims, self.dims)), dtype=self._matrix.dtype + self._backend.np.reshape(self._matrix, (self.dims, self.dims)), + dtype=self._matrix.dtype, ) reshaped = self._backend.cast( - np.transpose(np.conj(reshaped)) - reshaped, dtype=reshaped.dtype + self._backend.np.conj(reshaped).T - reshaped, dtype=reshaped.dtype ) norm = self._backend.calculate_norm_density_matrix(reshaped, order=order) @@ -240,19 +241,21 @@ def is_positive_semidefinite(self, precision_tol: float = 1e-8): self._matrix = self._full() self._pure = False - reshaped = np.reshape(self._matrix, (self.dims, self.dims)) - + reshaped = self._backend.np.reshape(self._matrix, (self.dims, self.dims)) if self.is_hermitian(): - eigenvalues = np.linalg.eigvalsh(reshaped) + eigenvalues = self._backend.calculate_eigenvalues(reshaped) else: if self._backend.__class__.__name__ in [ "CupyBackend", "CuQuantumBackend", ]: # pragma: no cover reshaped = np.array(reshaped.tolist(), dtype=reshaped.dtype) - eigenvalues = np.linalg.eigvals(reshaped) + eigenvalues = self._backend.calculate_eigenvalues(reshaped, hermitian=False) - return all(eigenvalue >= -precision_tol for eigenvalue in eigenvalues) + return all( + self._backend.np.real(eigenvalue) >= -precision_tol + for eigenvalue in eigenvalues + ) def is_channel( self, @@ -295,7 +298,9 @@ def apply(self, state): matrix = self._backend.cast(self._matrix, copy=True) if self.is_pure(): - return self._einsum("kj,ml,jl -> km", matrix, np.conj(matrix), state) + return self._einsum( + "kj,ml,jl -> km", matrix, self._backend.np.conj(matrix), state + ) return self._einsum("jklm,km -> jl", matrix, state) @@ -355,6 +360,7 @@ def link_product(self, second_network, subscripts: str = "ij,jk -> ik"): return QuantumNetwork( self._einsum(cexpr, first_matrix, second_matrix), [self.partition[0] + self.partition[-1]], + backend=self._backend, ) cexpr = "jkab,klbc->jlac" @@ -363,17 +369,19 @@ def link_product(self, second_network, subscripts: str = "ij,jk -> ik"): return QuantumNetwork( self._einsum(cexpr, second_matrix, first_matrix), [second_network.partition[0], self.partition[1]], + backend=self._backend, ) return QuantumNetwork( self._einsum(cexpr, first_matrix, second_matrix), [self.partition[0], second_network.partition[1]], + backend=self._backend, ) def copy(self): """Returns a copy of the :class:`qibo.quantum_info.quantum_networks.QuantumNetwork` object.""" return self.__class__( - np.copy(self._matrix), + self._backend.np.copy(self._matrix), partition=self.partition, system_output=self.system_output, pure=self._pure, @@ -640,10 +648,10 @@ def _set_tensor_and_parameters(self): try: if self._pure: - self._matrix = np.reshape(self._matrix, self.partition) + self._matrix = self._backend.np.reshape(self._matrix, self.partition) else: matrix_partition = self.partition * 2 - self._matrix = np.reshape(self._matrix, matrix_partition) + self._matrix = self._backend.np.reshape(self._matrix, matrix_partition) except: raise_error( ValueError, @@ -666,7 +674,9 @@ def _full(self): matrix = self._backend.cast(self._matrix, copy=True) if self.is_pure(): - matrix = self._einsum("jk,lm -> kjml", matrix, np.conj(matrix)) + matrix = self._einsum( + "jk,lm -> kjml", matrix, self._backend.np.conj(matrix) + ) return matrix diff --git a/src/qibo/quantum_info/random_ensembles.py b/src/qibo/quantum_info/random_ensembles.py index e6d9faf49b..59efb49b79 100644 --- a/src/qibo/quantum_info/random_ensembles.py +++ b/src/qibo/quantum_info/random_ensembles.py @@ -184,12 +184,12 @@ def random_hermitian( matrix = random_gaussian_matrix(dims, dims, seed=local_state, backend=backend) if semidefinite: - matrix = np.dot(np.transpose(np.conj(matrix)), matrix) + matrix = backend.np.matmul(backend.np.conj(matrix).T, matrix) else: - matrix = (matrix + np.transpose(np.conj(matrix))) / 2 + matrix = (matrix + backend.np.conj(matrix).T) / 2 if normalize: - matrix = matrix / np.linalg.norm(matrix) + matrix = matrix / np.linalg.norm(backend.to_numpy(matrix)) return matrix @@ -230,18 +230,18 @@ def random_unitary(dims: int, measure: Optional[str] = None, seed=None, backend= if measure == "haar": unitary = random_gaussian_matrix(dims, dims, seed=local_state, backend=backend) - Q, R = np.linalg.qr(unitary) - D = np.diag(R) - D = D / np.abs(D) - R = np.diag(D) - unitary = np.dot(Q, R) + # Tensorflow experi + Q, R = backend.np.linalg.qr(unitary) + D = backend.np.diag(R) + D = D / backend.np.abs(D) + R = backend.np.diag(D) + unitary = backend.np.matmul(Q, R) elif measure is None: from scipy.linalg import expm H = random_hermitian(dims, seed=seed, backend=NumpyBackend()) unitary = expm(-1.0j * H / 2) - - unitary = backend.cast(unitary, dtype=unitary.dtype) + unitary = backend.cast(unitary, dtype=unitary.dtype) return unitary @@ -351,7 +351,7 @@ def random_quantum_channel( else: super_op = random_unitary(dims, measure, local_state, backend) super_op = vectorization(super_op, order=order, backend=backend) - super_op = np.outer(super_op, np.conj(super_op)) + super_op = backend.np.outer(super_op, backend.np.conj(super_op)) if "chi" in representation: pauli_order = "IXYZ" @@ -436,10 +436,9 @@ def random_statevector(dims: int, seed=None, backend=None): backend, local_state = _check_backend_and_local_state(seed, backend) - state = local_state.standard_normal(dims).astype(complex) - state += 1.0j * local_state.standard_normal(dims) - state /= np.linalg.norm(state) - state = backend.cast(state, dtype=state.dtype) + state = backend.cast(local_state.standard_normal(dims).astype(complex)) + state = state + 1.0j * backend.cast(local_state.standard_normal(dims)) + state = state / backend.np.linalg.norm(state) return state @@ -546,24 +545,28 @@ def random_density_matrix( if pure: state = random_statevector(dims, seed=local_state, backend=backend) - state = np.outer(state, np.transpose(np.conj(state))) + state = backend.np.outer(state, backend.np.conj(state).T) else: if metric in ["hilbert-schmidt", "ginibre"]: state = random_gaussian_matrix( dims, rank, mean=0, stddev=1, seed=local_state, backend=backend ) - state = np.dot(state, np.transpose(np.conj(state))) - state = state / np.trace(state) + state = backend.np.matmul( + state, backend.np.transpose(backend.np.conj(state), (1, 0)) + ) + state = state / backend.np.trace(state) else: nqubits = int(np.log2(dims)) state = backend.identity_density_matrix(nqubits, normalize=False) state += random_unitary(dims, seed=local_state, backend=backend) - state = np.dot( + state = backend.np.matmul( state, random_gaussian_matrix(dims, rank, seed=local_state, backend=backend), ) - state = np.dot(state, np.transpose(np.conj(state))) - state /= np.trace(state) + state = backend.np.matmul( + state, backend.np.transpose(backend.np.conj(state), (1, 0)) + ) + state /= backend.np.trace(state) state = backend.cast(state, dtype=state.dtype) @@ -918,12 +921,15 @@ def random_pauli_hamiltonian( hamiltonian = random_hermitian(d, normalize=True, seed=local_state, backend=backend) - eigenvalues, eigenvectors = np.linalg.eigh(hamiltonian) + eigenvalues, eigenvectors = backend.calculate_eigenvectors(hamiltonian) + if backend.name == "tensorflow": + eigenvalues = backend.to_numpy(eigenvalues) + eigenvectors = backend.to_numpy(eigenvectors) if normalize is True: - eigenvalues -= eigenvalues[0] + eigenvalues = eigenvalues - eigenvalues[0] - eigenvalues /= eigenvalues[1] + eigenvalues = eigenvalues / eigenvalues[1] shift = 2 eigenvectors[:, shift:] = ( @@ -935,15 +941,17 @@ def random_pauli_hamiltonian( hamiltonian = backend.cast(hamiltonian, dtype=hamiltonian.dtype) # excluding the first eigenvector because first eigenvalue is zero for eigenvalue, eigenvector in zip( - eigenvalues[1:], np.transpose(eigenvectors)[1:] + eigenvalues[1:], backend.np.transpose(eigenvectors, (1, 0))[1:] ): - hamiltonian += eigenvalue * np.outer(eigenvector, np.conj(eigenvector)) + hamiltonian = hamiltonian + eigenvalue * backend.np.outer( + eigenvector, backend.np.conj(eigenvector) + ) U = comp_basis_to_pauli( nqubits, normalize=True, pauli_order=pauli_order, backend=backend ) - hamiltonian = np.real(U @ vectorization(hamiltonian, backend=backend)) + hamiltonian = backend.np.real(U @ vectorization(hamiltonian, backend=backend)) return hamiltonian, eigenvalues @@ -1193,10 +1201,12 @@ def _super_op_from_bcsz_measure(dims: int, rank: int, order: str, seed, backend) super_op = random_gaussian_matrix( dims**2, rank=rank, mean=0, stddev=1, seed=seed, backend=backend ) - super_op = super_op @ np.transpose(np.conj(super_op)) + super_op = super_op @ backend.np.conj(super_op).T # partial trace implemented with einsum - super_op_reduced = np.einsum("ijik->jk", np.reshape(super_op, (dims,) * 4)) + super_op_reduced = np.einsum( + "ijik->jk", np.reshape(backend.to_numpy(super_op), (dims,) * 4) + ) eigenvalues, eigenvectors = np.linalg.eigh(super_op_reduced) @@ -1204,8 +1214,12 @@ def _super_op_from_bcsz_measure(dims: int, rank: int, order: str, seed, backend) operator = np.zeros((dims, dims), dtype=complex) operator = backend.cast(operator, dtype=operator.dtype) - for eigenvalue, eigenvector in zip(eigenvalues, np.transpose(eigenvectors)): - operator += eigenvalue * np.outer(eigenvector, np.conj(eigenvector)) + for eigenvalue, eigenvector in zip( + backend.cast(eigenvalues), backend.cast(eigenvectors).T + ): + operator = operator + eigenvalue * backend.np.outer( + eigenvector, backend.np.conj(eigenvector) + ) if order == "row": operator = backend.np.kron( diff --git a/src/qibo/quantum_info/superoperator_transformations.py b/src/qibo/quantum_info/superoperator_transformations.py index 05c8622428..de1b76f2ec 100644 --- a/src/qibo/quantum_info/superoperator_transformations.py +++ b/src/qibo/quantum_info/superoperator_transformations.py @@ -64,12 +64,13 @@ def vectorization(state, order: str = "row", backend=None): backend = _check_backend(backend) if len(state.shape) == 1: - state = np.outer(state, np.conj(state)) + state = backend.np.outer(state, backend.np.conj(state)) if order == "row": - state = np.reshape(state, (1, -1), order="C")[0] + state = backend.np.reshape(state, (1, -1))[0] elif order == "column": - state = np.reshape(state, (1, -1), order="F")[0] + state = state.T + state = backend.np.reshape(state, (1, -1))[0] else: dim = len(state) nqubits = int(np.log2(dim)) @@ -78,11 +79,9 @@ def vectorization(state, order: str = "row", backend=None): for qubit in range(nqubits): new_axis += [qubit + nqubits, qubit] - state = np.reshape(state, [2] * 2 * nqubits) - state = np.transpose(state, axes=new_axis) - state = np.reshape(state, -1) - - state = backend.cast(state, dtype=state.dtype) + state = backend.np.reshape(state, [2] * 2 * nqubits) + state = backend.np.transpose(state, new_axis) + state = backend.np.reshape(state, (-1,)) return state @@ -130,20 +129,21 @@ def unvectorization(state, order: str = "row", backend=None): ) backend = _check_backend(backend) + state = backend.cast(state) dim = int(np.sqrt(len(state))) if order in ["row", "column"]: order = "C" if order == "row" else "F" - state = np.reshape(state, (dim, dim), order=order) + state = backend.cast( + np.reshape(backend.to_numpy(state), (dim, dim), order=order) + ) else: nqubits = int(np.log2(dim)) axes_old = list(np.arange(0, 2 * nqubits)) - state = np.reshape(state, [2] * 2 * nqubits) - state = np.transpose(state, axes=axes_old[1::2] + axes_old[0::2]) - state = np.reshape(state, [2**nqubits] * 2) - - state = backend.cast(state, dtype=state.dtype) + state = backend.np.reshape(state, [2] * 2 * nqubits) + state = backend.np.transpose(state, axes_old[1::2] + axes_old[0::2]) + state = backend.np.reshape(state, [2**nqubits] * 2) return state @@ -171,7 +171,7 @@ def to_choi(channel, order: str = "row", backend=None): ndarray: quantum channel in its Choi representation. """ channel = vectorization(channel, order=order, backend=backend) - channel = np.outer(channel, np.conj(channel)) + channel = backend.np.outer(channel, backend.np.conj(channel)) return channel @@ -239,7 +239,7 @@ def to_pauli_liouville( nqubits, normalize, pauli_order=pauli_order, backend=backend ) - channel = unitary @ channel @ np.transpose(np.conj(unitary)) + channel = unitary @ channel @ backend.np.conj(unitary).T return channel @@ -448,11 +448,12 @@ def choi_to_kraus( ) backend = _check_backend(backend) + choi_super_op = backend.cast(choi_super_op) if validate_cp: norm = float( backend.calculate_norm_density_matrix( - choi_super_op - np.transpose(np.conj(choi_super_op)), order=2 + choi_super_op - backend.np.conj(choi_super_op).T, order=2 ) ) if norm > PRECISION_TOL: @@ -460,22 +461,22 @@ def choi_to_kraus( else: # using eigh because, in this case, choi_super_op is # *already confirmed* to be Hermitian - eigenvalues, eigenvectors = np.linalg.eigh(choi_super_op) - eigenvectors = np.transpose(eigenvectors) + eigenvalues, eigenvectors = backend.calculate_eigenvectors(choi_super_op) + eigenvectors = eigenvectors.T - non_cp = bool(any(eigenvalues < -PRECISION_TOL)) + non_cp = bool(any(backend.np.real(eigenvalues) < -PRECISION_TOL)) else: non_cp = False # using eigh because, in this case, choi_super_op is # *assumed* to be Hermitian - eigenvalues, eigenvectors = np.linalg.eigh(choi_super_op) - eigenvectors = np.transpose(eigenvectors) + eigenvalues, eigenvectors = backend.calculate_eigenvectors(choi_super_op) + eigenvectors = eigenvectors.T if non_cp: warnings.warn("Input choi_super_op is a non-completely positive map.") # using singular value decomposition because choi_super_op is non-CP - U, coefficients, V = np.linalg.svd(choi_super_op) + U, coefficients, V = np.linalg.svd(backend.to_numpy(choi_super_op)) U = np.transpose(U) coefficients = np.sqrt(coefficients) V = np.conj(V) @@ -495,8 +496,8 @@ def choi_to_kraus( # when choi_super_op is CP kraus_ops, coefficients = [], [] for eig, kraus in zip(eigenvalues, eigenvectors): - if np.abs(eig) > precision_tol: - eig = np.sqrt(eig) + if backend.np.abs(eig) > precision_tol: + eig = backend.np.sqrt(eig) kraus_ops.append( eig * unvectorization(kraus, order=order, backend=backend) ) @@ -664,7 +665,7 @@ def kraus_to_choi(kraus_ops, order: str = "row", backend=None): kraus_op.append(gate) kraus_op = kraus_op.matrix(backend) kraus_op = vectorization(kraus_op, order=order, backend=backend) - super_op += np.outer(kraus_op, np.conj(kraus_op)) + super_op = super_op + backend.np.outer(kraus_op, backend.np.conj(kraus_op)) del kraus_op return super_op @@ -802,7 +803,7 @@ def kraus_to_chi( kraus_op = kraus_op.matrix(backend) kraus_op = vectorization(kraus_op, order=order, backend=backend) kraus_op = comp_to_pauli @ kraus_op - super_op += np.outer(kraus_op, np.conj(kraus_op)) + super_op = super_op + backend.np.outer(kraus_op, backend.np.conj(kraus_op)) del kraus_op return super_op @@ -872,7 +873,7 @@ def kraus_to_stinespring( # only utility is for outer product, # so np.conj here to only do it once - initial_state_env = np.conj(initial_state_env) + initial_state_env = backend.np.conj(initial_state_env) stinespring = np.zeros((dim_stinespring, dim_stinespring), dtype=complex) stinespring = backend.cast(stinespring, dtype=stinespring.dtype) @@ -884,9 +885,9 @@ def kraus_to_stinespring( kraus_op.append(gate) kraus_op = kraus_op.matrix(backend) kraus_op = backend.cast(kraus_op, dtype=kraus_op.dtype) - stinespring += np.kron( + stinespring = stinespring + backend.np.kron( kraus_op, - np.outer(vector_alpha, initial_state_env), + backend.np.outer(vector_alpha, initial_state_env), ) del kraus_op, vector_alpha @@ -979,7 +980,7 @@ def liouville_to_pauli( backend=backend, ) - return comp_to_pauli @ super_op @ np.conj(np.transpose(comp_to_pauli)) + return comp_to_pauli @ super_op @ backend.np.conj(comp_to_pauli.T) def liouville_to_kraus( @@ -1180,7 +1181,7 @@ def pauli_to_liouville( backend=backend, ) - return pauli_to_comp @ pauli_op @ np.conj(np.transpose(pauli_to_comp)) + return pauli_to_comp @ pauli_op @ backend.np.conj(pauli_to_comp).T def pauli_to_choi( @@ -1908,15 +1909,15 @@ def stinespring_to_kraus( initial_state_env, dtype=initial_state_env.dtype ) - stinespring = np.reshape(stinespring, (dim, dim_env, dim, dim_env)) - stinespring = np.swapaxes(stinespring, 1, 2) + stinespring = backend.np.reshape(stinespring, (dim, dim_env, dim, dim_env)) + stinespring = backend.np.swapaxes(stinespring, 1, 2) kraus_ops = [] for alpha in range(dim_env): vector_alpha = np.zeros(dim_env, dtype=complex) vector_alpha[alpha] = 1.0 vector_alpha = backend.cast(vector_alpha, dtype=vector_alpha.dtype) - kraus = np.conj(vector_alpha) @ stinespring @ initial_state_env + kraus = backend.np.conj(vector_alpha) @ stinespring @ initial_state_env kraus_ops.append(kraus) return kraus_ops @@ -2057,7 +2058,7 @@ def function(x0, operators): operator = (1 - np.sum(x0)) * np.eye(dim**2, dtype=complex) operator = backend.cast(operator, dtype=operator.dtype) for prob, oper in zip(x0, operators): - operator += prob * oper + operator = operator + prob * oper return float(backend.calculate_norm_density_matrix(target - operator, order=2)) @@ -2123,6 +2124,7 @@ def _reshuffling(super_op, order: str = "row", backend=None): Returns: ndarray: Choi (Liouville) representation of the quantum channel. """ + super_op = backend.cast(super_op) if not isinstance(order, str): raise_error(TypeError, f"order must be type str, but it is type {type(order)}.") @@ -2152,13 +2154,12 @@ def _reshuffling(super_op, order: str = "row", backend=None): raise_error(ValueError, "super_op must be of shape (4^n, 4^n)") dim = int(dim) - super_op = np.reshape(super_op, [dim] * 4) + super_op = backend.np.reshape(super_op, [dim] * 4) axes = [1, 2] if order == "row" else [0, 3] - super_op = np.swapaxes(super_op, *axes) + super_op = backend.np.swapaxes(super_op, *axes) - super_op = np.reshape(super_op, [dim**2, dim**2]) - super_op = backend.cast(super_op, dtype=super_op.dtype) + super_op = backend.np.reshape(super_op, [dim**2, dim**2]) return super_op @@ -2217,7 +2218,7 @@ def _individual_kraus_to_liouville( kraus_op.append(gate) kraus_op = kraus_op.matrix(backend) kraus_op = vectorization(kraus_op, order=order, backend=backend) - kraus_op = np.outer(kraus_op, np.conj(kraus_op)) + kraus_op = backend.np.outer(kraus_op, backend.np.conj(kraus_op)) super_ops.append(choi_to_liouville(kraus_op, order=order, backend=backend)) return super_ops diff --git a/src/qibo/quantum_info/utils.py b/src/qibo/quantum_info/utils.py index 38ae21f012..ef21f2c97f 100644 --- a/src/qibo/quantum_info/utils.py +++ b/src/qibo/quantum_info/utils.py @@ -185,11 +185,10 @@ def hadamard_transform(array, implementation: str = "fast", backend=None): return array - array = _hadamard_transform_1d(array) + array = _hadamard_transform_1d(array, backend=backend) if len(array.shape) == 2: - array = _hadamard_transform_1d(np.transpose(array)) - array = np.transpose(array) + array = _hadamard_transform_1d(array.T, backend=backend).T # needed for the tensorflow backend array = backend.cast(array, dtype=array.dtype) @@ -252,7 +251,10 @@ def hellinger_distance(prob_dist_p, prob_dist_q, validate: bool = False, backend raise_error(ValueError, "Second probability array must sum to 1.") distance = float( - backend.calculate_norm(np.sqrt(prob_dist_p) - np.sqrt(prob_dist_q)) / np.sqrt(2) + backend.calculate_norm( + backend.np.sqrt(prob_dist_p) - backend.np.sqrt(prob_dist_q) + ) + / np.sqrt(2) ) return distance @@ -396,11 +398,15 @@ def haar_integral( rand_unit_density, dtype=rand_unit_density.dtype ) for _ in range(samples): - haar_state = np.reshape(random_statevector(dim, backend=backend), (-1, 1)) + haar_state = backend.np.reshape( + random_statevector(dim, backend=backend), (-1, 1) + ) - rho = haar_state @ np.conj(np.transpose(haar_state)) + rho = haar_state @ backend.np.conj(haar_state).T - rand_unit_density += reduce(np.kron, [rho] * power_t) + rand_unit_density = rand_unit_density + reduce( + backend.np.kron, [rho] * power_t + ) integral = rand_unit_density / samples @@ -414,14 +420,16 @@ def haar_integral( ] identity = np.eye(dim**power_t, dtype=float) - identity = backend.cast(identity, dtype=identity.dtype) identity = np.reshape(identity, (dim,) * (2 * power_t)) + identity = backend.cast(identity, dtype=identity.dtype) integral = np.zeros((dim**power_t, dim**power_t), dtype=float) integral = backend.cast(integral, dtype=integral.dtype) for indices in permutations_list: - integral += np.reshape(np.transpose(identity, indices), (-1, dim**power_t)) - integral *= normalization + integral = integral + backend.np.reshape( + backend.np.transpose(identity, indices), (-1, dim**power_t) + ) + integral = integral * normalization return integral @@ -469,25 +477,26 @@ def pqc_integral(circuit, power_t: int, samples: int, backend=None): rho = backend.execute_circuit(circuit).state() - rand_unit_density += reduce(np.kron, [rho] * power_t) + rand_unit_density = rand_unit_density + reduce(np.kron, [rho] * power_t) integral = rand_unit_density / samples return integral -def _hadamard_transform_1d(array): +def _hadamard_transform_1d(array, backend=None): # necessary because of tf.EagerTensor # does not accept item assignment - array_copied = np.copy(array) + backend = _check_backend(backend) + array_copied = backend.np.copy(array) indexes = [2**k for k in range(int(np.log2(len(array_copied))))] for index in indexes: for k in range(0, len(array_copied), 2 * index): for j in range(k, k + index): # copy necessary because of cupy backend - elem_1 = np.copy(array_copied[j]) - elem_2 = np.copy(array_copied[j + index]) + elem_1 = backend.np.copy(array_copied[j]) + elem_2 = backend.np.copy(array_copied[j + index]) array_copied[j] = elem_1 + elem_2 array_copied[j + index] = elem_1 - elem_2 array_copied /= 2.0 diff --git a/src/qibo/result.py b/src/qibo/result.py index a10ea75076..905cbb8e49 100644 --- a/src/qibo/result.py +++ b/src/qibo/result.py @@ -288,7 +288,9 @@ def probabilities(self, qubits: Optional[Union[list, set]] = None): probs[state] = freq / self.nshots probs = self.backend.cast(probs) self._probs = probs - return self.backend.calculate_probabilities(np.sqrt(probs), qubits, nqubits) + return self.backend.calculate_probabilities( + self.backend.np.sqrt(probs), qubits, nqubits + ) def has_samples(self): """Check whether the samples are available already. @@ -323,7 +325,11 @@ def samples(self, binary: bool = True, registers: bool = False): if self._samples is None: if self.measurements[0].result.has_samples(): self._samples = self.backend.np.concatenate( - [gate.result.samples() for gate in self.measurements], axis=1 + [ + self.backend.cast(gate.result.samples()) + for gate in self.measurements + ], + axis=1, ) else: if self._frequencies is not None: @@ -350,7 +356,7 @@ def samples(self, binary: bool = True, registers: bool = False): qubit_map = { q: i for i, q in enumerate(self.measurement_gate.target_qubits) } - self._samples = self.backend.cast(samples, "int32") + self._samples = samples for gate in self.measurements: rqubits = tuple(qubit_map.get(q) for q in gate.target_qubits) gate.result.register_samples( diff --git a/src/qibo/transpiler/unitary_decompositions.py b/src/qibo/transpiler/unitary_decompositions.py index cf8509b623..ad892db734 100644 --- a/src/qibo/transpiler/unitary_decompositions.py +++ b/src/qibo/transpiler/unitary_decompositions.py @@ -1,8 +1,8 @@ import numpy as np from qibo import gates, matrices -from qibo.backends import NumpyBackend, _check_backend -from qibo.config import raise_error +from qibo.backends import _check_backend +from qibo.config import PRECISION_TOL, raise_error magic_basis = np.array( [[1, -1j, 0, 0], [0, 0, 1, -1j], [0, 0, -1, -1j], [1, 1j, 0, 0]] @@ -48,7 +48,7 @@ def calculate_psi(unitary, magic_basis=magic_basis, backend=None): backend (:class:`qibo.backends.abstract.Backend`): Backend to use for calculations. Returns: - (ndarray) Eigenvectors in the computational basis and eigenvalues of :math:`U^{T} U`. + ndarray: Eigenvectors in the computational basis and eigenvalues of :math:`U^{T} U`. """ backend = _check_backend(backend) @@ -71,15 +71,14 @@ def calculate_psi(unitary, magic_basis=magic_basis, backend=None): ) # construct and diagonalize UT_U ut_u = backend.np.transpose(u_magic, (1, 0)) @ u_magic - # When the matrix given to np.linalg.eig is a diagonal matrix up to machine precision the decomposition - # is not accurate anymore. decimals = 20 works for random 2q Clifford unitaries. - if backend.__class__.__name__ == "TensorflowBackend": + if backend.__class__.__name__ != "PyTorchBackend": + # 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)) - psi_magic, _ = np.linalg.qr(psi_magic) else: - eigvals, psi_magic = backend.np.linalg.eig(np.round(ut_u, decimals=20)) - # orthogonalize eigenvectors in the case of degeneracy (Gram-Schmidt) - psi_magic, _ = backend.np.linalg.qr(psi_magic) + eigvals, psi_magic = backend.calculate_eigenvectors(ut_u, hermitian=False) + # orthogonalize eigenvectors in the case of degeneracy (Gram-Schmidt) + psi_magic, _ = backend.np.linalg.qr(psi_magic) # write psi in computational basis psi = backend.np.matmul(magic_basis, psi_magic) return psi, eigvals @@ -96,11 +95,12 @@ def schmidt_decompose(state, backend=None): """ backend = _check_backend(backend) + # tf.linalg.svd has a different behaviour if backend.__class__.__name__ == "TensorflowBackend": u, d, v = np.linalg.svd(backend.np.reshape(state, (2, 2))) else: u, d, v = backend.np.linalg.svd(backend.np.reshape(state, (2, 2))) - if not np.allclose(d, [1, 0]): # pragma: no cover + if not np.allclose(backend.to_numpy(d), [1, 0]): # pragma: no cover raise_error( ValueError, f"Unexpected singular values: {d}\nCan only decompose product states.", @@ -121,8 +121,11 @@ def calculate_single_qubit_unitaries(psi, backend=None): """ backend = _check_backend(backend) psi_magic = backend.np.matmul(backend.np.conj(backend.cast(magic_basis)).T, psi) - if not np.allclose( - backend.to_numpy(psi_magic).imag, np.zeros_like(psi_magic) + if ( + backend.np.real( + backend.calculate_norm_density_matrix(backend.np.imag(psi_magic)) + ) + > PRECISION_TOL ): # pragma: no cover raise_error(NotImplementedError, "Given state is not real in the magic basis.") psi_bar = backend.cast(psi.T, copy=True) @@ -210,10 +213,10 @@ def to_bell_diagonal(ud, bell_basis=bell_basis, backend=None): backend.np.transpose(backend.np.conj(bell_basis), (1, 0)) @ ud @ bell_basis ) ud_diag = backend.np.diag(ud_bell) - if not np.allclose(backend.np.diag(ud_diag), ud_bell): # pragma: no cover + if not backend.np.allclose(backend.np.diag(ud_diag), ud_bell): # pragma: no cover return None uprod = backend.np.prod(ud_diag) - if not np.allclose(uprod, 1): # pragma: no cover + if not np.allclose(backend.to_numpy(uprod), 1): # pragma: no cover return None return ud_diag diff --git a/tests/test_backends.py b/tests/test_backends.py index aa250db686..41b5c99df3 100644 --- a/tests/test_backends.py +++ b/tests/test_backends.py @@ -128,3 +128,19 @@ def test_list_available_backends(): "qibotn": {"cutensornet": False, "qutensornet": True}, } assert available_backends == list_available_backends() + + +def test_gradients_pytorch(): + from qibo.backends import PyTorchBackend # pylint: disable=import-outside-toplevel + + backend = PyTorchBackend() + gate = gates.RX(0, 0.1) + matrix = gate.matrix(backend) + assert matrix.requires_grad + assert backend.gradients + backend.requires_grad(False) + gate = gates.RX(0, 0.1) + matrix = gate.matrix(backend) + assert not matrix.requires_grad + assert not backend.gradients + assert not backend.matrices.requires_grad diff --git a/tests/test_callbacks.py b/tests/test_callbacks.py index bf2ac7d2d0..a551245b1a 100644 --- a/tests/test_callbacks.py +++ b/tests/test_callbacks.py @@ -98,7 +98,9 @@ def test_entropy_in_circuit(backend, density_matrix, base): backend.assert_allclose(values, target, atol=PRECISION_TOL) target_spectrum = [0.0] + list([0, 0, np.log(2), np.log(2)] / np.log(base)) - entropy_spectrum = np.ravel(np.concatenate(entropy.spectrum)).tolist() + entropy_spectrum = backend.np.ravel( + backend.np.concatenate(entropy.spectrum) + ).tolist() backend.assert_allclose(entropy_spectrum, target_spectrum, atol=PRECISION_TOL) @@ -170,7 +172,7 @@ def target_entropy(t): target = [0, target_entropy(0.1234), 0, target_entropy(0.4321)] values = [backend.to_numpy(x) for x in entropy[:]] - backend.assert_allclose(values, target, atol=PRECISION_TOL) + backend.assert_allclose(values, target, atol=1e-6) c = Circuit(8, accelerators) with pytest.raises(RuntimeError): @@ -228,17 +230,21 @@ def test_entropy_density_matrix(backend): rho = random_density_matrix(2**4, backend=backend) # this rho is not always positive. Make rho positive for this application - _, u = np.linalg.eigh(rho) + _, u = np.linalg.eigh(backend.to_numpy(rho)) + u = backend.cast(u, dtype=u.dtype) matrix = np.random.random(u.shape[0]) - matrix = backend.cast(matrix, dtype=matrix.dtype) - rho = np.dot(np.dot(u, np.diag(5 * matrix)), np.conj(np.transpose(u))) + matrix = backend.cast(matrix, dtype=u.dtype) + rho = backend.np.matmul( + backend.np.matmul(u, backend.np.diag(5 * matrix)), + backend.np.conj(backend.np.transpose(u, (1, 0))), + ) # this is a positive rho entropy = callbacks.EntanglementEntropy([1, 3]) entropy.nqubits = 4 final_ent = entropy.apply(backend, rho) - rho = rho.reshape(8 * (2,)) + rho = backend.to_numpy(rho.reshape(8 * (2,))) reduced_rho = np.einsum("abcdafch->bdfh", rho).reshape((4, 4)) eigvals = np.linalg.eigvalsh(reduced_rho).real # assert that all eigenvalues are non-negative @@ -278,12 +284,12 @@ def test_norm(backend, density_matrix, seed): if density_matrix: norm.nqubits = 1 state = random_density_matrix(2**norm.nqubits, seed=seed, backend=backend) - target_norm = np.trace(state) + target_norm = backend.np.trace(state) final_norm = norm.apply_density_matrix(backend, state) else: norm.nqubits = 2 state = random_statevector(2**norm.nqubits, seed=seed, backend=backend) - target_norm = np.sqrt((np.abs(state) ** 2).sum()) + target_norm = np.sqrt((np.abs(backend.to_numpy(state)) ** 2).sum()) final_norm = norm.apply(backend, state) backend.assert_allclose(final_norm, target_norm) @@ -330,13 +336,16 @@ def test_energy(backend, density_matrix): from qibo.quantum_info import random_density_matrix state = random_density_matrix(2**4, backend=backend) - target_energy = np.trace(np.dot(matrix, state)) + target_energy = backend.np.trace(backend.np.matmul(matrix, state)) final_energy = energy.apply_density_matrix(backend, state) else: from qibo.quantum_info import random_statevector state = random_statevector(2**4, backend=backend) - target_energy = np.dot(np.conj(state), np.dot(matrix, state)) + target_energy = np.matmul( + np.conj(backend.to_numpy(state)), + np.matmul(backend.to_numpy(matrix), backend.to_numpy(state)), + ) final_energy = energy.apply(backend, state) backend.assert_allclose(final_energy, target_energy) @@ -356,7 +365,7 @@ def test_gap(backend, dense, check_degenerate): ham = lambda t: (1 - t) * h0.matrix + t * h1.matrix targets = {"ground": [], "excited": [], "gap": []} for t in np.linspace(0, 1, 11): - eigvals = np.real(np.linalg.eigvalsh(ham(t))) + eigvals = np.real(np.linalg.eigvalsh(backend.to_numpy(ham(t)))) targets["ground"].append(eigvals[0]) targets["excited"].append(eigvals[1]) targets["gap"].append(eigvals[1] - eigvals[0]) diff --git a/tests/test_derivative.py b/tests/test_derivative.py index a1c6e97e0e..8a39efcee7 100644 --- a/tests/test_derivative.py +++ b/tests/test_derivative.py @@ -2,7 +2,6 @@ import pytest from qibo import Circuit, gates, hamiltonians -from qibo.backends.pytorch import PyTorchBackend from qibo.derivative import finite_differences, parameter_shift from qibo.symbols import Z @@ -47,55 +46,31 @@ def test_standard_parameter_shift(backend, nshots, atol, scale_factor, grads): # testing parameter out of bounds with pytest.raises(ValueError): - grad_0 = parameter_shift( + grad = parameter_shift( circuit=c, hamiltonian=test_hamiltonian, parameter_index=5 ) # testing hamiltonian type with pytest.raises(TypeError): - grad_0 = parameter_shift( + grad = parameter_shift( circuit=c, hamiltonian=c, parameter_index=0, nshots=nshots ) - if isinstance(backend, PyTorchBackend): - with pytest.raises(NotImplementedError) as excinfo: - grad = parameter_shift( - circuit=c, hamiltonian=test_hamiltonian, parameter_index=0 - ) - assert ( - str(excinfo.value) - == "PyTorchBackend for the parameter shift rule is not supported." - ) - - else: - # executing all the procedure - grad_0 = parameter_shift( - circuit=c, - hamiltonian=test_hamiltonian, - parameter_index=0, - scale_factor=scale_factor, - nshots=nshots, - ) - grad_1 = parameter_shift( - circuit=c, - hamiltonian=test_hamiltonian, - parameter_index=1, - scale_factor=scale_factor, - nshots=nshots, - ) - grad_2 = parameter_shift( + # executing all the procedure + grad = [ + parameter_shift( circuit=c, hamiltonian=test_hamiltonian, - parameter_index=2, + parameter_index=i, scale_factor=scale_factor, nshots=nshots, ) + for i in range(3) + ] - # check of known values - # calculated using tf.GradientTape - backend.assert_allclose(grad_0, grads[0], atol=atol) - backend.assert_allclose(grad_1, grads[1], atol=atol) - backend.assert_allclose(grad_2, grads[2], atol=atol) + # check of known values + for i in range(3): + backend.assert_allclose(grad[i], grads[i], atol=atol) @pytest.mark.parametrize("step_size", [10**-i for i in range(5, 10, 1)]) diff --git a/tests/test_gates_channels.py b/tests/test_gates_channels.py index a512b83831..def0cd9a30 100644 --- a/tests/test_gates_channels.py +++ b/tests/test_gates_channels.py @@ -14,21 +14,28 @@ def test_general_channel(backend): """""" - a_1 = np.sqrt(0.4) * matrices.X - a_2 = np.sqrt(0.6) * np.array( - [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]] + a_1 = backend.cast(np.sqrt(0.4) * matrices.X) + a_2 = backend.cast( + np.sqrt(0.6) + * np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) ) initial_state = random_density_matrix(2**2, backend=backend) - m_1 = np.kron(np.eye(2), backend.to_numpy(a_1)) + m_1 = backend.np.kron(backend.identity_density_matrix(1, normalize=False), a_1) m_1 = backend.cast(m_1, dtype=m_1.dtype) m_2 = backend.cast(a_2, dtype=a_2.dtype) - target_state = np.dot(np.dot(m_1, initial_state), np.transpose(np.conj(m_1))) - target_state += np.dot(np.dot(m_2, initial_state), np.transpose(np.conj(m_2))) + target_state = backend.np.matmul( + backend.np.matmul(m_1, initial_state), + backend.np.transpose(backend.np.conj(m_1), (1, 0)), + ) + target_state = target_state + backend.np.matmul( + backend.np.matmul(m_2, initial_state), + backend.np.transpose(backend.np.conj(m_2), (1, 0)), + ) channel1 = gates.KrausChannel([(1,), (0, 1)], [a_1, a_2]) assert channel1.target_qubits == (0, 1) final_state = backend.apply_channel_density_matrix( - channel1, np.copy(initial_state), 2 + channel1, backend.np.copy(initial_state), 2 ) backend.assert_allclose(final_state, target_state) @@ -37,7 +44,7 @@ def test_general_channel(backend): channel2 = gates.KrausChannel([(1,), (0, 1)], [a_1, a_2]) assert channel2.target_qubits == (0, 1) final_state = backend.apply_channel_density_matrix( - channel2, np.copy(initial_state), 2 + channel2, backend.np.copy(initial_state), 2 ) backend.assert_allclose(final_state, target_state) @@ -85,10 +92,12 @@ def test_kraus_channel(backend, pauli_order): [0.4 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.6 + 0.0j], ] ) - test_choi = np.reshape(test_superop, [2] * 4).swapaxes(0, 3).reshape([4, 4]) + test_choi = backend.cast( + np.reshape(test_superop, [2] * 4).swapaxes(0, 3).reshape([4, 4]) + ) pauli_elements = {"I": 2.0, "X": -0.4, "Y": -2.0, "Z": 0.4} - test_pauli = np.diag([pauli_elements[p] for p in pauli_order]) + test_pauli = backend.cast(np.diag([pauli_elements[p] for p in pauli_order])) test_superop = backend.cast(test_superop, dtype=test_superop.dtype) test_choi = backend.cast(test_choi, dtype=test_choi.dtype) @@ -130,8 +139,10 @@ def test_kraus_channel(backend, pauli_order): def test_unitary_channel(backend): """""" - a_1 = matrices.X - a_2 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) + a_1 = backend.cast(matrices.X) + a_2 = backend.cast( + np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) + ) qubits = [(0,), (2, 3)] probabilities = [0.4, 0.3] @@ -140,18 +151,18 @@ def test_unitary_channel(backend): channel = gates.UnitaryChannel(qubits, matrices_) initial_state = random_density_matrix(2**4, backend=backend) final_state = backend.apply_channel_density_matrix( - channel, np.copy(initial_state), 4 + channel, backend.np.copy(initial_state), 4 ) - eye = np.eye(2) - ma_1 = np.kron(np.kron(a_1, eye), np.kron(eye, eye)) - ma_2 = np.kron(np.kron(eye, eye), a_2) + eye = backend.identity_density_matrix(1, normalize=False) + ma_1 = backend.np.kron(backend.np.kron(a_1, eye), backend.np.kron(eye, eye)) + ma_2 = backend.np.kron(backend.np.kron(eye, eye), a_2) ma_1 = backend.cast(ma_1, dtype=ma_1.dtype) ma_2 = backend.cast(ma_2, dtype=ma_2.dtype) target_state = ( 0.3 * initial_state - + 0.4 * np.dot(ma_1, np.dot(initial_state, ma_1)) - + 0.3 * np.dot(ma_2, np.dot(initial_state, ma_2)) + + 0.4 * backend.np.matmul(ma_1, backend.np.matmul(initial_state, ma_1)) + + 0.3 * backend.np.matmul(ma_2, backend.np.matmul(initial_state, ma_2)) ) backend.assert_allclose(final_state, target_state) @@ -174,7 +185,7 @@ def test_unitary_channel_probability_tolerance(backend): matrices_ = [(p, np.random.random((4, 4))) for p in probs] identity_channel = gates.UnitaryChannel(qubits, matrices_) backend.assert_allclose( - identity_channel.to_liouville(backend=backend), np.eye(num_terms) + identity_channel.to_liouville(backend=backend), backend.np.eye(num_terms) ) @@ -205,10 +216,12 @@ def test_pauli_noise_channel(backend, pauli_order): qubits = (1,) channel = gates.PauliNoiseChannel(qubits, [("X", 0.3)]) final_state = backend.apply_channel_density_matrix( - channel, np.copy(initial_state), 2 + channel, backend.np.copy(initial_state), 2 ) gate = gates.X(1) - target_state = backend.apply_gate_density_matrix(gate, np.copy(initial_state), 2) + target_state = backend.apply_gate_density_matrix( + gate, backend.np.copy(initial_state), 2 + ) target_state = 0.3 * target_state + 0.7 * initial_state backend.assert_allclose(final_state, target_state) @@ -248,7 +261,9 @@ def test_depolarizing_channel(backend): lam = 0.3 initial_state_r = backend.partial_trace_density_matrix(initial_state, (2,), 3) channel = gates.DepolarizingChannel((0, 1), lam) - final_state = channel.apply_density_matrix(backend, np.copy(initial_state), 3) + final_state = channel.apply_density_matrix( + backend, backend.np.copy(initial_state), 3 + ) final_state_r = backend.partial_trace_density_matrix(final_state, (2,), 3) target_state_r = (1 - lam) * initial_state_r + lam * backend.cast( np.identity(4) @@ -272,10 +287,12 @@ def test_amplitude_damping_channel(backend): channel = gates.AmplitudeDampingChannel(0, gamma) initial_state = random_density_matrix(2**1, backend=backend) - final_state = channel.apply_density_matrix(backend, np.copy(initial_state), 1) - target_state = kraus_0 @ initial_state @ np.transpose( - np.conj(kraus_0) - ) + kraus_1 @ initial_state @ np.transpose(np.conj(kraus_1)) + final_state = channel.apply_density_matrix( + backend, backend.np.copy(initial_state), 1 + ) + target_state = kraus_0 @ initial_state @ backend.np.transpose( + backend.np.conj(kraus_0), (1, 0) + ) + kraus_1 @ initial_state @ backend.np.transpose(backend.np.conj(kraus_1), (1, 0)) backend.assert_allclose(final_state, target_state) @@ -296,10 +313,12 @@ def test_phase_damping_channel(backend): channel = gates.PhaseDampingChannel(0, gamma) initial_state = random_density_matrix(2**1, backend=backend) - final_state = channel.apply_density_matrix(backend, np.copy(initial_state), 1) - target_state = kraus_0 @ initial_state @ np.transpose( - np.conj(kraus_0) - ) + kraus_1 @ initial_state @ np.transpose(np.conj(kraus_1)) + final_state = channel.apply_density_matrix( + backend, backend.np.copy(initial_state), 1 + ) + target_state = kraus_0 @ initial_state @ backend.np.transpose( + backend.np.conj(kraus_0), (1, 0) + ) + kraus_1 @ initial_state @ backend.np.transpose(backend.np.conj(kraus_1), (1, 0)) backend.assert_allclose(final_state, target_state) @@ -311,7 +330,7 @@ def test_thermal_relaxation_channel(backend, t_1, t_2, time, excpop): """Check ``gates.ThermalRelaxationChannel`` on a 3-qubit random density matrix.""" initial_state = random_density_matrix(2**3, backend=backend) gate = gates.ThermalRelaxationChannel(0, [t_1, t_2, time, excpop]) - final_state = gate.apply_density_matrix(backend, np.copy(initial_state), 3) + final_state = gate.apply_density_matrix(backend, backend.np.copy(initial_state), 3) if t_2 > t_1: p_0, p_1, exp = ( @@ -323,8 +342,10 @@ def test_thermal_relaxation_channel(backend, t_1, t_2, time, excpop): matrix[0, -1], matrix[-1, 0] = exp, exp matrix = matrix.reshape(4 * (2,)) # Apply matrix using Eq. (3.28) from arXiv:1111.6950 - target_state = np.copy(initial_state).reshape(6 * (2,)) - target_state = np.einsum("abcd,aJKcjk->bJKdjk", matrix, target_state) + target_state = backend.np.copy(initial_state).reshape(6 * (2,)) + target_state = np.einsum( + "abcd,aJKcjk->bJKdjk", matrix, backend.to_numpy(target_state) + ) target_state = target_state.reshape(initial_state.shape) else: p_0, p_1, p_z = ( @@ -332,7 +353,10 @@ def test_thermal_relaxation_channel(backend, t_1, t_2, time, excpop): gate.init_kwargs["p_1"], gate.init_kwargs["p_z"], ) - m_z = np.kron(matrices.Z, np.kron(matrices.I, matrices.I)) + m_z = backend.np.kron( + backend.cast(matrices.Z), + backend.np.kron(backend.cast(matrices.I), backend.cast(matrices.I)), + ) m_z = backend.cast(m_z, dtype=m_z.dtype) z_rho = m_z @ initial_state @ m_z @@ -353,7 +377,9 @@ def test_thermal_relaxation_channel(backend, t_1, t_2, time, excpop): ones = backend.cast(ones, dtype=ones.dtype) target_state = (1 - p_0 - p_1 - p_z) * initial_state + p_z * z_rho - target_state += np.reshape(p_0 * zeros + p_1 * ones, initial_state.shape) + target_state += backend.np.reshape( + p_0 * zeros + p_1 * ones, initial_state.shape + ) target_state = backend.cast(target_state, dtype=target_state.dtype) @@ -398,7 +424,7 @@ def test_readout_error_channel(backend): probability_sum = gates.ReadoutErrorChannel( 0, stochastic_noise ).apply_density_matrix(backend, rho, 1) - probability_sum = np.diag(probability_sum).sum().real + probability_sum = np.diag(backend.to_numpy(probability_sum)).sum().real backend.assert_allclose(probability_sum - 1 < PRECISION_TOL, True) @@ -407,7 +433,9 @@ def test_reset_channel(backend): """""" initial_state = random_density_matrix(2**3, backend=backend) gate = gates.ResetChannel(0, [0.2, 0.2]) - final_state = backend.reset_error_density_matrix(gate, np.copy(initial_state), 3) + final_state = backend.reset_error_density_matrix( + gate, backend.np.copy(initial_state), 3 + ) trace = backend.to_numpy( backend.partial_trace_density_matrix(initial_state, (0,), 3) @@ -421,7 +449,7 @@ def test_reset_channel(backend): zeros = backend.cast(zeros, dtype=zeros.dtype) ones = backend.cast(ones, dtype=ones.dtype) - target_state = 0.6 * initial_state + 0.2 * np.reshape( + target_state = 0.6 * initial_state + 0.2 * backend.np.reshape( zeros + ones, initial_state.shape ) diff --git a/tests/test_gates_density_matrix.py b/tests/test_gates_density_matrix.py index b83e02a98e..ec4fde3e27 100644 --- a/tests/test_gates_density_matrix.py +++ b/tests/test_gates_density_matrix.py @@ -9,7 +9,7 @@ def apply_gates(backend, gatelist, nqubits=None, initial_state=None): - state = backend.cast(np.copy(initial_state)) + state = backend.cast(backend.np.copy(initial_state)) for gate in gatelist: state = backend.apply_gate_density_matrix(gate, state, nqubits) return backend.to_numpy(state) @@ -21,10 +21,10 @@ def test_hgate_density_matrix(backend): gate = gates.H(1) final_rho = apply_gates(backend, [gate], 2, initial_rho) - matrix = np.array([[1, 1], [1, -1]]) / np.sqrt(2) + matrix = np.array([[1 + 0j, 1], [1 + 0j, -1]]) / np.sqrt(2) matrix = np.kron(np.eye(2), matrix) matrix = backend.cast(matrix, dtype=matrix.dtype) - target_rho = np.dot(np.dot(matrix, initial_rho), matrix) + target_rho = backend.np.matmul(backend.np.matmul(matrix, initial_rho), matrix) backend.assert_allclose(final_rho, target_rho) @@ -38,7 +38,9 @@ def test_rygate_density_matrix(backend): phase = np.exp(1j * theta / 2.0) matrix = phase * np.array([[phase.real, -phase.imag], [phase.imag, phase.real]]) matrix = backend.cast(matrix, dtype=matrix.dtype) - target_rho = np.dot(np.dot(matrix, initial_rho), np.transpose(np.conj(matrix))) + target_rho = backend.np.matmul( + backend.np.matmul(matrix, initial_rho), backend.np.conj(matrix).T + ) backend.assert_allclose(final_rho, target_rho, atol=PRECISION_TOL) @@ -72,7 +74,9 @@ def test_one_qubit_gates(backend, gatename, gatekwargs): final_rho = apply_gates(backend, [gate], 1, initial_rho) matrix = backend.to_numpy(gate.matrix(backend)) - target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T) + target_rho = np.einsum( + "ab,bc,cd->ad", matrix, backend.to_numpy(initial_rho), matrix.conj().T + ) backend.assert_allclose(final_rho, target_rho) @@ -81,14 +85,16 @@ def test_controlled_by_one_qubit_gates(backend, gatename): nqubits = 2 initial_rho = random_density_matrix(2**nqubits, seed=1, backend=backend) gate = getattr(gates, gatename)(1).controlled_by(0) - final_rho = apply_gates(backend, [gate], 2, np.copy(initial_rho)) + final_rho = apply_gates(backend, [gate], 2, backend.np.copy(initial_rho)) matrix = backend.to_numpy(backend.matrix(getattr(gates, gatename)(1))) cmatrix = np.eye(4, dtype=matrix.dtype) cmatrix[2:, 2:] = matrix - cmatrix = backend.cast(cmatrix, dtype=cmatrix.dtype) target_rho = np.einsum( - "ab,bc,cd->ad", cmatrix, initial_rho, np.transpose(np.conj(cmatrix)) + "ab,bc,cd->ad", + cmatrix, + backend.to_numpy(initial_rho), + np.transpose(np.conj(cmatrix)), ) backend.assert_allclose(final_rho, target_rho) @@ -116,9 +122,9 @@ def test_two_qubit_gates(backend, gatename, gatekwargs): gate = getattr(gates, gatename)(0, 1, **gatekwargs) final_rho = apply_gates(backend, [gate], 2, initial_rho) - matrix = gate.matrix(backend) + matrix = backend.to_numpy(gate.matrix(backend)) target_rho = np.einsum( - "ab,bc,cd->ad", matrix, initial_rho, np.transpose(np.conj(matrix)) + "ab,bc,cd->ad", matrix, backend.to_numpy(initial_rho), np.conj(matrix).T ) backend.assert_allclose(final_rho, target_rho, atol=PRECISION_TOL) @@ -131,7 +137,9 @@ def test_toffoli_gate(backend): final_rho = apply_gates(backend, [gate], 3, initial_rho) matrix = backend.to_numpy(gate.matrix(backend)) - target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T) + target_rho = np.einsum( + "ab,bc,cd->ad", matrix, backend.to_numpy(initial_rho), np.conj(matrix).T + ) backend.assert_allclose(final_rho, target_rho) @@ -143,7 +151,9 @@ def test_unitary_gate(backend, nqubits): initial_rho = random_density_matrix(2**nqubits, backend=backend) gate = gates.Unitary(matrix, *range(nqubits)) final_rho = apply_gates(backend, [gate], nqubits, initial_rho) - target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T) + target_rho = np.einsum( + "ab,bc,cd->ad", matrix, backend.to_numpy(initial_rho), np.conj(matrix).T + ) backend.assert_allclose(final_rho, target_rho) @@ -159,7 +169,9 @@ def test_cu1gate_application_twoqubit(backend): matrix[3, 3] = np.exp(1j * theta) matrix = np.kron(matrix, np.eye(2)) matrix = backend.cast(matrix, dtype=matrix.dtype) - target_rho = np.dot(np.dot(matrix, initial_rho), np.transpose(np.conj(matrix))) + target_rho = backend.np.matmul( + backend.np.matmul(matrix, initial_rho), backend.np.conj(matrix).T + ) backend.assert_allclose(final_rho, target_rho) @@ -190,13 +202,17 @@ def test_controlled_with_effect(backend): c.add(gates.X(0)) c.add(gates.X(2)) c.add(gates.SWAP(1, 3).controlled_by(0, 2)) - final_rho = backend.execute_circuit(c, np.copy(initial_rho)).state() + final_rho = backend.execute_circuit( + c, backend.np.copy(backend.cast(initial_rho)) + ).state() c = Circuit(4, density_matrix=True) c.add(gates.X(0)) c.add(gates.X(2)) c.add(gates.SWAP(1, 3)) - target_rho = backend.execute_circuit(c, np.copy(initial_rho)).state() + target_rho = backend.execute_circuit( + c, backend.np.copy(backend.cast(initial_rho)) + ).state() backend.assert_allclose(final_rho, target_rho) @@ -205,15 +221,15 @@ def test_controlled_by_random(backend, nqubits): """Check controlled_by method on gate.""" initial_psi = random_statevector(2**nqubits, backend=backend) - initial_rho = np.outer(initial_psi, np.conj(initial_psi)) + initial_rho = backend.np.outer(initial_psi, backend.np.conj(initial_psi)) c = Circuit(nqubits, density_matrix=True) c.add(gates.RX(1, theta=0.789).controlled_by(2)) c.add(gates.fSim(0, 2, theta=0.123, phi=0.321).controlled_by(1, 3)) - final_rho = backend.execute_circuit(c, np.copy(initial_rho)).state() + final_rho = backend.execute_circuit(c, backend.np.copy(initial_rho)).state() c = Circuit(nqubits) c.add(gates.RX(1, theta=0.789).controlled_by(2)) c.add(gates.fSim(0, 2, theta=0.123, phi=0.321).controlled_by(1, 3)) - target_psi = backend.execute_circuit(c, np.copy(initial_psi)).state() - target_rho = np.outer(target_psi, np.conj(target_psi)) + target_psi = backend.execute_circuit(c, backend.np.copy(initial_psi)).state() + target_rho = backend.np.outer(target_psi, backend.np.conj(target_psi)) backend.assert_allclose(final_rho, target_rho) diff --git a/tests/test_gates_gates.py b/tests/test_gates_gates.py index c0a3cd10ec..27028e4a4c 100644 --- a/tests/test_gates_gates.py +++ b/tests/test_gates_gates.py @@ -27,8 +27,8 @@ def apply_gates(backend, gatelist, nqubits=None, initial_state=None): def test_h(backend): final_state = apply_gates(backend, [gates.H(0), gates.H(1)], nqubits=2) - target_state = np.ones_like(final_state) / 2 - backend.assert_allclose(final_state, target_state) + target_state = np.ones_like(backend.to_numpy(final_state)) / 2 + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.H(0).qasm_label == "h" assert gates.H(0).clifford @@ -37,9 +37,9 @@ def test_h(backend): def test_x(backend): final_state = apply_gates(backend, [gates.X(0)], nqubits=2) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[2] = 1.0 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.X(0).qasm_label == "x" assert gates.X(0).clifford @@ -48,9 +48,9 @@ def test_x(backend): def test_y(backend): final_state = apply_gates(backend, [gates.Y(1)], nqubits=2) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[1] = 1j - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.Y(0).qasm_label == "y" assert gates.Y(0).clifford @@ -59,10 +59,10 @@ def test_y(backend): def test_z(backend): final_state = apply_gates(backend, [gates.H(0), gates.H(1), gates.Z(0)], nqubits=2) - target_state = np.ones_like(final_state) / 2.0 + target_state = np.ones_like(backend.to_numpy(final_state)) / 2.0 target_state[2] *= -1.0 target_state[3] *= -1.0 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.Z(0).qasm_label == "z" assert gates.Z(0).clifford @@ -86,11 +86,11 @@ def test_sx(backend): initial_state=initial_state, ) - matrix = np.array([[1 + 1j, 1 - 1j], [1 - 1j, 1 + 1j]]) / 2 - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = np.array([[1 + 1j, 1 - 1j], [1 - 1j, 1 + 1j]]) / 2.0 + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference observable = random_hermitian(2**nqubits, backend=backend) @@ -98,10 +98,9 @@ def test_sx(backend): np_obs = backend.to_numpy(observable) np_target_state = backend.to_numpy(target_state) backend.assert_allclose( - np.transpose(np.conj(np_final_state_decompose)) - @ np_obs - @ np_final_state_decompose, - np.transpose(np.conj(np_target_state)) @ np_obs @ np_target_state, + np.conj(np_final_state_decompose).T @ np_obs @ np_final_state_decompose, + np.conj(np_target_state).T @ np_obs @ np_target_state, + atol=1e-6, ) assert gates.SX(0).qasm_label == "sx" @@ -127,10 +126,10 @@ def test_sxdg(backend): ) matrix = np.array([[1 - 1j, 1 + 1j], [1 + 1j, 1 - 1j]]) / 2 - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference observable = random_hermitian(2**nqubits, backend=backend) @@ -138,10 +137,9 @@ def test_sxdg(backend): np_obs = backend.to_numpy(observable) np_target_state = backend.to_numpy(target_state) backend.assert_allclose( - np.transpose(np.conj(np_final_state_decompose)) - @ np_obs - @ np_final_state_decompose, - np.transpose(np.conj(np_target_state)) @ np_obs @ np_target_state, + np.conj(np_final_state_decompose).T @ np_obs @ np_final_state_decompose, + np.conj(np_target_state).T @ np_obs @ np_target_state, + atol=1e-6, ) assert gates.SXDG(0).qasm_label == "sxdg" @@ -152,7 +150,7 @@ def test_sxdg(backend): def test_s(backend): final_state = apply_gates(backend, [gates.H(0), gates.H(1), gates.S(1)], nqubits=2) target_state = np.array([0.5, 0.5j, 0.5, 0.5j]) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.S(0).qasm_label == "s" assert gates.S(0).clifford @@ -164,7 +162,7 @@ def test_sdg(backend): backend, [gates.H(0), gates.H(1), gates.SDG(1)], nqubits=2 ) target_state = np.array([0.5, -0.5j, 0.5, -0.5j]) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.SDG(0).qasm_label == "sdg" assert gates.SDG(0).clifford @@ -174,7 +172,7 @@ def test_sdg(backend): def test_t(backend): final_state = apply_gates(backend, [gates.H(0), gates.H(1), gates.T(1)], nqubits=2) target_state = np.array([0.5, (1 + 1j) / np.sqrt(8), 0.5, (1 + 1j) / np.sqrt(8)]) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.T(0).qasm_label == "t" assert not gates.T(0).clifford @@ -186,7 +184,7 @@ def test_tdg(backend): backend, [gates.H(0), gates.H(1), gates.TDG(1)], nqubits=2 ) target_state = np.array([0.5, (1 - 1j) / np.sqrt(8), 0.5, (1 - 1j) / np.sqrt(8)]) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.TDG(0).qasm_label == "tdg" assert not gates.TDG(0).clifford @@ -196,11 +194,11 @@ def test_tdg(backend): def test_identity(backend): gatelist = [gates.H(0), gates.H(1), gates.I(0), gates.I(1)] final_state = apply_gates(backend, gatelist, nqubits=2) - target_state = np.ones_like(final_state) / 2.0 - backend.assert_allclose(final_state, target_state) + target_state = np.ones_like(backend.to_numpy(final_state)) / 2.0 + backend.assert_allclose(final_state, target_state, atol=1e-6) gatelist = [gates.H(0), gates.H(1), gates.I(0, 1)] final_state = apply_gates(backend, gatelist, nqubits=2) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.I(0).qasm_label == "id" assert gates.I(0).clifford @@ -221,11 +219,11 @@ def test_align(backend): final_state = apply_gates(backend, gate_list, nqubits=nqubits) target_state = backend.plus_state(nqubits) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) gate_matrix = gate.matrix(backend) identity = backend.identity_density_matrix(nqubits, normalize=False) - backend.assert_allclose(gate_matrix, identity) + backend.assert_allclose(gate_matrix, identity, atol=1e-6) with pytest.raises(NotImplementedError): gate.qasm_label @@ -237,7 +235,7 @@ def test_align(backend): # :class:`qibo.core.cgates.M` is tested seperately in `test_measurement_gate.py` -@pytest.mark.parametrize("theta", [np.random.rand(), np.pi / 2, -np.pi / 2, np.pi]) +@pytest.mark.parametrize("theta", [np.random.rand(), np.pi / 2.0, -np.pi / 2.0, np.pi]) def test_rx(backend, theta): nqubits = 1 initial_state = random_statevector(2**nqubits, backend=backend) @@ -250,14 +248,14 @@ def test_rx(backend, theta): phase = np.exp(1j * theta / 2.0) gate = np.array([[phase.real, -1j * phase.imag], [-1j * phase.imag, phase.real]]) - gate = backend.cast(gate, dtype=gate.dtype) + gate = backend.cast(gate) target_state = gate @ initial_state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.RX(0, theta=theta).qasm_label == "rx" assert gates.RX(0, theta=theta).unitary - if (theta % (np.pi / 2)).is_integer(): + if (theta % (np.pi / 2.0)).is_integer(): assert gates.RX(0, theta=theta).clifford else: assert not gates.RX(0, theta=theta).clifford @@ -274,7 +272,7 @@ def test_rx(backend, theta): ) -@pytest.mark.parametrize("theta", [np.random.rand(), np.pi / 2, -np.pi / 2, np.pi]) +@pytest.mark.parametrize("theta", [np.random.rand(), np.pi / 2.0, -np.pi / 2.0, np.pi]) def test_ry(backend, theta): nqubits = 1 initial_state = random_statevector(2**nqubits, backend=backend) @@ -287,21 +285,21 @@ def test_ry(backend, theta): phase = np.exp(1j * theta / 2.0) gate = np.array([[phase.real, -phase.imag], [phase.imag, phase.real]]) - gate = backend.cast(gate, dtype="complex128") - target_state = gate @ backend.cast(initial_state, dtype="complex128") + gate = backend.cast(gate) + target_state = gate @ backend.cast(initial_state) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.RY(0, theta=theta).qasm_label == "ry" assert gates.RY(0, theta=theta).unitary - if (theta % (np.pi / 2)).is_integer(): + if (theta % (np.pi / 2.0)).is_integer(): assert gates.RY(0, theta=theta).clifford else: assert not gates.RY(0, theta=theta).clifford @pytest.mark.parametrize("apply_x", [True, False]) -@pytest.mark.parametrize("theta", [np.random.rand(), np.pi / 2, -np.pi / 2, np.pi]) +@pytest.mark.parametrize("theta", [np.random.rand(), np.pi / 2.0, -np.pi / 2.0, np.pi]) def test_rz(backend, theta, apply_x): nqubits = 1 if apply_x: @@ -310,10 +308,10 @@ def test_rz(backend, theta, apply_x): gatelist = [] gatelist.append(gates.RZ(0, theta)) final_state = apply_gates(backend, gatelist, nqubits=nqubits) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) p = int(apply_x) target_state[p] = np.exp((2 * p - 1) * 1j * theta / 2.0) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.RZ(0, theta).qasm_label == "rz" assert gates.RZ(0, theta=theta).unitary @@ -367,10 +365,10 @@ def test_gpi(backend): phase = np.exp(1.0j * phi) matrix = np.array([[0, np.conj(phase)], [phase, 0]]) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) - target_state = np.dot(matrix, initial_state) - backend.assert_allclose(final_state, target_state) + target_state = backend.np.matmul(matrix, initial_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.GPI(0, phi).qasm_label == "gpi" assert not gates.GPI(0, phi).clifford @@ -387,10 +385,10 @@ def test_gpi2(backend): phase = np.exp(1.0j * phi) matrix = np.array([[1, -1.0j * np.conj(phase)], [-1.0j * phase, 1]]) / np.sqrt(2) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) - target_state = np.dot(matrix, initial_state) - backend.assert_allclose(final_state, target_state) + target_state = backend.np.matmul(matrix, initial_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.GPI2(0, phi).qasm_label == "gpi2" assert not gates.GPI2(0, phi).clifford @@ -400,9 +398,9 @@ def test_gpi2(backend): def test_u1(backend): theta = 0.1234 final_state = apply_gates(backend, [gates.X(0), gates.U1(0, theta)], nqubits=1) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[1] = np.exp(1j * theta) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.U1(0, theta).qasm_label == "u1" assert not gates.U1(0, theta).clifford @@ -423,11 +421,9 @@ def test_u2(backend): [np.exp(1j * (phi - lam) / 2), np.exp(1j * (phi + lam) / 2)], ] ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + target_state = np.matmul(matrix, backend.to_numpy(initial_state)) / np.sqrt(2) - target_state = np.dot(matrix, initial_state) / np.sqrt(2) - - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.U2(0, phi, lam).qasm_label == "u2" assert not gates.U2(0, phi, lam).clifford @@ -457,21 +453,22 @@ def test_u3(backend, seed_state, seed_observable): em = np.exp(1j * (phi - lam) / 2) matrix = np.array([[ep.conj() * cost, -em.conj() * sint], [em * sint, ep * cost]]) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) - target_state = np.dot(matrix, initial_state) + target_state = backend.np.matmul(matrix, initial_state) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference observable = random_hermitian(2**nqubits, seed=seed_observable, backend=backend) backend.assert_allclose( - backend.cast(np.transpose(np.conj(final_state_decompose))) + backend.cast(backend.np.conj(final_state_decompose).T) @ observable @ final_state_decompose, - backend.cast(np.transpose(np.conj(target_state))) + backend.cast(backend.np.conj(target_state).T) @ observable @ backend.cast(target_state), + atol=1e-6, ) assert gates.U3(0, theta, phi, lam).qasm_label == "u3" assert not gates.U3(0, theta, phi, lam).clifford @@ -492,11 +489,11 @@ def test_u1q(backend, seed_state): matrix = np.array( [[cost, -1j * np.exp(-1j * phi) * sint], [-1j * np.exp(1j * phi) * sint, cost]] ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) - target_state = np.dot(matrix, initial_state) + target_state = backend.np.matmul(matrix, initial_state) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert not gates.U1q(0, theta, phi).clifford assert gates.U1q(0, theta, phi).unitary @@ -510,9 +507,9 @@ def test_cnot(backend, applyx): gatelist = [] gatelist.append(gates.CNOT(0, 1)) final_state = apply_gates(backend, gatelist, nqubits=2) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[3 * int(applyx)] = 1.0 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.CNOT(0, 1).qasm_label == "cx" assert gates.CNOT(0, 1).clifford @@ -532,9 +529,9 @@ def test_cy(backend, seed_state, seed_observable): [0, 0, 1j, 0], ] ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) - target_state = np.dot(matrix, initial_state) + target_state = backend.np.matmul(matrix, initial_state) # test decomposition final_state_decompose = apply_gates( backend, @@ -547,17 +544,18 @@ def test_cy(backend, seed_state, seed_observable): final_state = apply_gates(backend, [gate], initial_state=initial_state) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference observable = random_hermitian(2**nqubits, seed=seed_observable, backend=backend) backend.assert_allclose( - backend.cast(np.transpose(np.conj(final_state_decompose))) + backend.cast(backend.np.conj(final_state_decompose).T) @ observable @ final_state_decompose, - backend.cast(np.transpose(np.conj(target_state))) + backend.cast(backend.np.conj(target_state).T) @ observable @ backend.cast(target_state), + atol=1e-6, ) assert gate.name == "cy" @@ -573,9 +571,9 @@ def test_cz(backend, seed_state, seed_observable): initial_state = random_statevector(2**nqubits, seed=seed_state, backend=backend) matrix = np.eye(4) matrix[3, 3] = -1 - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) - target_state = np.dot(matrix, initial_state) + target_state = backend.np.matmul(matrix, initial_state) # test decomposition final_state_decompose = apply_gates( backend, @@ -588,17 +586,18 @@ def test_cz(backend, seed_state, seed_observable): final_state = apply_gates(backend, [gate], initial_state=initial_state) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference observable = random_hermitian(2**nqubits, seed=seed_observable, backend=backend) backend.assert_allclose( - backend.cast(np.transpose(np.conj(final_state_decompose))) + backend.cast(backend.np.conj(final_state_decompose).T) @ observable @ final_state_decompose, - backend.cast(np.transpose(np.conj(target_state))) + backend.cast(backend.np.conj(target_state).T) @ observable @ backend.cast(target_state), + atol=1e-6, ) assert gate.name == "cz" @@ -632,11 +631,11 @@ def test_csx(backend): [0, 0, (1 - 1j) / 2, (1 + 1j) / 2], ] ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) - backend.assert_allclose(final_state_decompose, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + backend.assert_allclose(final_state_decompose, target_state, atol=1e-6) assert gates.CSX(0, 1).qasm_label == "csx" assert not gates.CSX(0, 1).clifford @@ -668,11 +667,11 @@ def test_csxdg(backend): [0, 0, (1 + 1j) / 2, (1 - 1j) / 2], ] ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) - backend.assert_allclose(final_state_decompose, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + backend.assert_allclose(final_state_decompose, target_state, atol=1e-6) assert gates.CSXDG(0, 1).qasm_label == "csxdg" assert not gates.CSXDG(0, 1).clifford @@ -719,16 +718,16 @@ def test_cun(backend, name, params): _matrix = gate.matrix(backend) gate = backend.cast(_matrix, dtype=_matrix.dtype) - target_state = np.dot(gate, initial_state) + target_state = backend.np.matmul(gate, initial_state) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_swap(backend): final_state = apply_gates(backend, [gates.X(1), gates.SWAP(0, 1)], nqubits=2) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[2] = 1.0 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.SWAP(0, 1).qasm_label == "swap" assert gates.SWAP(0, 1).clifford @@ -737,9 +736,9 @@ def test_swap(backend): def test_iswap(backend): final_state = apply_gates(backend, [gates.X(1), gates.iSWAP(0, 1)], nqubits=2) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[2] = 1.0j - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.iSWAP(0, 1).qasm_label == "iswap" assert gates.iSWAP(0, 1).clifford @@ -748,10 +747,10 @@ def test_iswap(backend): def test_siswap(backend): final_state = apply_gates(backend, [gates.X(1), gates.SiSWAP(0, 1)], nqubits=2) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[1] = 1.0 / np.sqrt(2) target_state[2] = 1.0j / np.sqrt(2) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert not gates.SiSWAP(0, 1).clifford assert gates.SiSWAP(0, 1).unitary @@ -783,11 +782,11 @@ def test_fswap(backend): ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) - backend.assert_allclose(final_state_decompose, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + backend.assert_allclose(final_state_decompose, target_state, atol=1e-6) assert gates.FSWAP(0, 1).qasm_label == "fswap" assert gates.FSWAP(0, 1).clifford @@ -799,7 +798,7 @@ def test_multiple_swap(backend): final_state = apply_gates(backend, gatelist, nqubits=4) gatelist = [gates.X(1), gates.X(3)] target_state = apply_gates(backend, gatelist, nqubits=4) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_fsim(backend): @@ -807,17 +806,17 @@ def test_fsim(backend): phi = 0.4321 gatelist = [gates.H(0), gates.H(1), gates.fSim(0, 1, theta, phi)] final_state = apply_gates(backend, gatelist, nqubits=2) - target_state = np.ones_like(final_state) / 2.0 + target_state = np.ones_like(backend.to_numpy(final_state)) / 2.0 rotation = np.array( [[np.cos(theta), -1j * np.sin(theta)], [-1j * np.sin(theta), np.cos(theta)]] ) matrix = np.eye(4, dtype=target_state.dtype) matrix[1:3, 1:3] = rotation matrix[3, 3] = np.exp(-1j * phi) - matrix = backend.cast(matrix, dtype=matrix.dtype) - target_state = np.dot(matrix, target_state) + matrix = backend.cast(matrix) + target_state = backend.np.matmul(matrix, backend.cast(target_state)) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gates.fSim(0, 1, theta, phi).qasm_label @@ -838,17 +837,17 @@ def test_sycamore(backend): matrix = np.array( [ - [1, 0, 0, 0], + [1 + 0j, 0, 0, 0], [0, 0, -1j, 0], [0, -1j, 0, 0], [0, 0, 0, np.exp(-1j * np.pi / 6)], ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gates.SYC(0, 1).qasm_label @@ -863,14 +862,15 @@ def test_generalized_fsim(backend): gatelist = [gates.H(0), gates.H(1), gates.H(2)] gatelist.append(gates.GeneralizedfSim(1, 2, rotation, phi)) final_state = apply_gates(backend, gatelist, nqubits=3) - target_state = np.ones_like(final_state) / np.sqrt(8) + target_state = np.ones(len(final_state), dtype=complex) / np.sqrt(8) matrix = np.eye(4, dtype=target_state.dtype) matrix[1:3, 1:3] = rotation matrix[3, 3] = np.exp(-1j * phi) - matrix = backend.cast(matrix, dtype=matrix.dtype) - target_state[:4] = np.dot(matrix, target_state[:4]) - target_state[4:] = np.dot(matrix, target_state[4:]) - backend.assert_allclose(final_state, target_state) + # matrix = backend.cast(matrix, dtype=matrix.dtype) + target_state[:4] = np.matmul(matrix, target_state[:4]) + target_state[4:] = np.matmul(matrix, target_state[4:]) + target_state = backend.cast(target_state, dtype=target_state.dtype) + backend.assert_allclose(final_state, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gatelist[-1].qasm_label @@ -883,7 +883,7 @@ def test_generalized_fsim_parameter_setter(backend): phi = np.random.random() matrix = np.random.random((2, 2)) gate = gates.GeneralizedfSim(0, 1, matrix, phi) - backend.assert_allclose(gate.parameters[0], matrix) + backend.assert_allclose(gate.parameters[0], matrix, atol=1e-6) assert gate.parameters[1] == phi @@ -914,7 +914,7 @@ def test_rxx(backend): ] ) target_state = gate.dot(np.ones(4)) / 2.0 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.RXX(0, 1, theta=theta).qasm_label == "rxx" assert not gates.RXX(0, 1, theta).clifford @@ -936,7 +936,7 @@ def test_ryy(backend): ] ) target_state = gate.dot(np.ones(4)) / 2.0 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.RYY(0, 1, theta=theta).qasm_label == "ryy" assert not gates.RYY(0, 1, theta).clifford @@ -948,9 +948,9 @@ def test_rzz(backend): final_state = apply_gates( backend, [gates.X(0), gates.X(1), gates.RZZ(0, 1, theta=theta)], nqubits=2 ) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[3] = np.exp(-1j * theta / 2.0) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.RZZ(0, 1, theta=theta).qasm_label == "rzz" assert not gates.RZZ(0, 1, theta).clifford @@ -985,11 +985,11 @@ def test_rzx(backend): ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) - backend.assert_allclose(final_state_decompose, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + backend.assert_allclose(final_state_decompose, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gates.RZX(0, 1, theta).qasm_label @@ -1026,18 +1026,19 @@ def test_rxxyy(backend): ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state observable = random_hermitian(2**nqubits, backend=backend) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference backend.assert_allclose( - backend.cast(np.transpose(np.conj(final_state_decompose))) + backend.cast(backend.np.conj(final_state_decompose).T) @ observable @ final_state_decompose, - backend.cast(np.transpose(np.conj(target_state))) @ observable @ target_state, + backend.cast(backend.np.conj(target_state).T) @ observable @ target_state, + atol=1e-6, ) with pytest.raises(NotImplementedError): @@ -1060,7 +1061,7 @@ def test_ms(backend): [gates.H(0), gates.H(1), gates.MS(0, 1, phi0=phi0, phi1=phi1)], nqubits=2, ) - target_state = np.ones_like(final_state) / 2.0 + target_state = np.ones_like(backend.to_numpy(final_state)) / 2.0 plus = np.exp(1.0j * (phi0 + phi1)) minus = np.exp(1.0j * (phi0 - phi1)) @@ -1070,10 +1071,10 @@ def test_ms(backend): matrix[2, 1] = -1.0j * minus matrix[1, 2] = -1.0j * np.conj(minus) matrix /= np.sqrt(2) - matrix = backend.cast(matrix, dtype=matrix.dtype) - target_state = np.dot(matrix, target_state) + matrix = backend.cast(matrix) + target_state = backend.np.matmul(matrix, backend.cast(target_state)) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gates.MS(0, 1, phi0, phi1, theta).qasm_label == "ms" assert not gates.MS(0, 1, phi0, phi1).clifford @@ -1107,11 +1108,11 @@ def test_givens(backend): ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) - backend.assert_allclose(final_state_decompose, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + backend.assert_allclose(final_state_decompose, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gates.GIVENS(0, 1, theta).qasm_label @@ -1147,11 +1148,11 @@ def test_rbs(backend): ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) - backend.assert_allclose(final_state_decompose, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + backend.assert_allclose(final_state_decompose, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gates.RBS(0, 1, theta).qasm_label @@ -1186,17 +1187,18 @@ def test_ecr(backend): ], dtype=np.complex128, ) / np.sqrt(2) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # testing random expectation value due to global phase difference observable = random_hermitian(2**nqubits, backend=backend) backend.assert_allclose( - backend.cast(np.transpose(np.conj(final_state_decompose))) + backend.cast(backend.np.conj(final_state_decompose).T) @ observable @ final_state_decompose, - backend.cast(np.transpose(np.conj(target_state))) @ observable @ target_state, + backend.cast(backend.np.conj(target_state).T) @ observable @ target_state, + atol=1e-6, ) with pytest.raises(NotImplementedError): @@ -1213,12 +1215,12 @@ def test_toffoli(backend, applyx): else: gatelist = [gates.X(1), gates.TOFFOLI(0, 1, 2)] final_state = apply_gates(backend, gatelist, nqubits=3) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) if applyx: target_state[-1] = 1 else: target_state[2] = 1 - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) assert gatelist[-1].qasm_label == "ccx" assert not gates.TOFFOLI(0, 1, 2).clifford @@ -1297,10 +1299,10 @@ def test_deutsch(backend): ], dtype=np.complex128, ) - matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = backend.cast(matrix) target_state = matrix @ initial_state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) with pytest.raises(NotImplementedError): gates.DEUTSCH(0, 1, 2, theta).qasm_label @@ -1318,14 +1320,14 @@ def test_unitary(backend, nqubits): gate = gates.Unitary(matrix, *range(1, nqubits), name="random") gatelist.append(gate) final_state = apply_gates(backend, gatelist, nqubits=nqubits) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_unitary_initialization(backend): matrix = np.random.random((4, 4)) gate = gates.Unitary(matrix, 0, 1) - backend.assert_allclose(gate.parameters[0], matrix) + backend.assert_allclose(gate.parameters[0], matrix, atol=1e-6) with pytest.raises(NotImplementedError): gates.Unitary(matrix, 0, 1).qasm_label @@ -1342,7 +1344,7 @@ def test_unitary_common_gates(backend): gates.Unitary(np.array([[1, 1], [1, -1]]) / np.sqrt(2), 1), ] final_state = apply_gates(backend, gatelist, nqubits=2) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) thetax = 0.1234 thetay = 0.4321 @@ -1364,7 +1366,7 @@ def test_unitary_common_gates(backend): cnot = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) gatelist = [gates.Unitary(rx, 0), gates.Unitary(ry, 1), gates.Unitary(cnot, 0, 1)] final_state = apply_gates(backend, gatelist, nqubits=2) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_unitary_multiqubit(backend): @@ -1382,7 +1384,7 @@ def test_unitary_multiqubit(backend): unitary = gates.Unitary(matrix, 0, 1, 2, 3) final_state = apply_gates(backend, [unitary], nqubits=4) target_state = apply_gates(backend, gatelist, nqubits=4) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) ############################# Test ``controlled_by`` ############################# @@ -1400,7 +1402,7 @@ def test_controlled_x(backend): final_state = apply_gates(backend, gatelist, nqubits=4) gatelist = [gates.X(1), gates.X(3)] target_state = apply_gates(backend, gatelist, nqubits=4) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_controlled_x_vs_cnot(backend): @@ -1408,7 +1410,7 @@ def test_controlled_x_vs_cnot(backend): final_state = apply_gates(backend, gatelist, nqubits=3) gatelist = [gates.X(0), gates.CNOT(0, 2)] target_state = apply_gates(backend, gatelist, nqubits=3) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_controlled_x_vs_toffoli(backend): @@ -1416,7 +1418,7 @@ def test_controlled_x_vs_toffoli(backend): final_state = apply_gates(backend, gatelist, nqubits=3) gatelist = [gates.X(0), gates.X(2), gates.TOFFOLI(0, 2, 1)] target_state = apply_gates(backend, gatelist, nqubits=3) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) @pytest.mark.parametrize("applyx", [False, True]) @@ -1434,7 +1436,7 @@ def test_controlled_rx(backend, applyx): gatelist.extend([gates.X(1), gates.RX(2, theta)]) target_state = apply_gates(backend, gatelist, nqubits=3) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_controlled_u1(backend): @@ -1444,9 +1446,11 @@ def test_controlled_u1(backend): gatelist.append(gates.X(0)) gatelist.append(gates.X(1)) final_state = apply_gates(backend, gatelist, nqubits=3) - target_state = np.zeros_like(final_state) + target_state = np.zeros_like(backend.to_numpy(final_state)) target_state[1] = np.exp(1j * theta) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) + gate = gates.U1(0, theta).controlled_by(1) + assert gate.__class__.__name__ == "CU1" def test_controlled_u2(backend): @@ -1458,7 +1462,7 @@ def test_controlled_u2(backend): final_state = apply_gates(backend, gatelist, nqubits=3) gatelist = [gates.X(0), gates.X(1), gates.U2(2, phi, lam), gates.X(0), gates.X(1)] target_state = apply_gates(backend, gatelist, nqubits=3) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # for coverage gate = gates.CU2(0, 1, phi, lam) assert gate.parameters == (phi, lam) @@ -1474,7 +1478,7 @@ def test_controlled_u3(backend): target_state = apply_gates( backend, [gates.CU3(0, 1, theta, phi, lam)], 2, initial_state ) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) # for coverage gate = gates.U3(0, theta, phi, lam) assert gate.parameters == (theta, phi, lam) @@ -1499,7 +1503,7 @@ def test_controlled_swap(backend, applyx, free_qubit): if applyx: gatelist.extend([gates.X(0), gates.SWAP(1 + f, 2 + f)]) target_state = apply_gates(backend, gatelist, 3 + f) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) @pytest.mark.parametrize("applyx", [False, True]) @@ -1517,7 +1521,7 @@ def test_controlled_swap_double(backend, applyx): if applyx: gatelist.extend([gates.X(3), gates.SWAP(1, 2)]) target_state = apply_gates(backend, gatelist, 4) - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_controlled_fsim(backend): @@ -1526,40 +1530,39 @@ def test_controlled_fsim(backend): gatelist.append(gates.fSim(5, 3, theta, phi).controlled_by(0, 2, 1)) final_state = apply_gates(backend, gatelist, 6) - target_state = np.ones_like(final_state) / np.sqrt(2**6) + target_state = np.ones(len(final_state), dtype=complex) / np.sqrt(2**6) rotation = np.array( [[np.cos(theta), -1j * np.sin(theta)], [-1j * np.sin(theta), np.cos(theta)]] ) matrix = np.eye(4, dtype=target_state.dtype) matrix[1:3, 1:3] = rotation matrix[3, 3] = np.exp(-1j * phi) - matrix = backend.cast(matrix, dtype=matrix.dtype) ids = [56, 57, 60, 61] - target_state[ids] = np.dot(matrix, target_state[ids]) + target_state[ids] = np.matmul(matrix, target_state[ids]) ids = [58, 59, 62, 63] - target_state[ids] = np.dot(matrix, target_state[ids]) - backend.assert_allclose(final_state, target_state) + target_state[ids] = np.matmul(matrix, target_state[ids]) + target_state = backend.cast(target_state, dtype=target_state.dtype) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_controlled_unitary(backend): - matrix = np.random.random((2, 2)) - # matrix = backend.cast(matrix, dtype=matrix.dtype) + matrix = random_unitary(2**1, backend=backend) gatelist = [gates.H(0), gates.H(1), gates.Unitary(matrix, 1).controlled_by(0)] final_state = apply_gates(backend, gatelist, 2) - target_state = np.ones_like(final_state) / 2.0 - matrix = backend.cast(matrix, dtype=matrix.dtype) - target_state[2:] = np.dot(matrix, target_state[2:]) - backend.assert_allclose(final_state, target_state) + target_state = np.ones(len(final_state), dtype=complex) / 2.0 + target_state[2:] = np.matmul(backend.to_numpy(matrix), target_state[2:]) + target_state = backend.cast(target_state, dtype=target_state.dtype) + backend.assert_allclose(final_state, target_state, atol=1e-6) - matrix = np.random.random((4, 4)) + matrix = random_unitary(2**2, backend=backend) gatelist = [gates.H(i) for i in range(4)] gatelist.append(gates.Unitary(matrix, 1, 3).controlled_by(0, 2)) final_state = apply_gates(backend, gatelist, 4) - target_state = np.ones_like(final_state) / 4.0 + target_state = np.ones(len(final_state), dtype=complex) / 4.0 ids = [10, 11, 14, 15] - matrix = backend.cast(matrix, dtype=matrix.dtype) - target_state[ids] = np.dot(matrix, target_state[ids]) - backend.assert_allclose(final_state, target_state) + target_state[ids] = np.matmul(backend.to_numpy(matrix), target_state[ids]) + target_state = backend.cast(target_state, dtype=target_state.dtype) + backend.assert_allclose(final_state, target_state, atol=1e-6) ############################################################################### @@ -1617,7 +1620,7 @@ def test_dagger(backend, gate, args): nqubits = len(gate.qubits) initial_state = random_statevector(2**nqubits, backend=backend) final_state = apply_gates(backend, [gate, gate.dagger()], nqubits, initial_state) - backend.assert_allclose(final_state, initial_state) + backend.assert_allclose(final_state, initial_state, atol=1e-6) GATES = [ @@ -1640,7 +1643,7 @@ def test_controlled_dagger(backend, gate, args): nqubits = 4 initial_state = random_statevector(2**nqubits, backend=backend) final_state = apply_gates(backend, [gate, gate.dagger()], 4, initial_state) - backend.assert_allclose(final_state, initial_state) + backend.assert_allclose(final_state, initial_state, atol=1e-6) @pytest.mark.parametrize("gate_1,gate_2", [("S", "SDG"), ("T", "TDG")]) @@ -1650,19 +1653,19 @@ def test_dagger_consistency(backend, gate_1, gate_2, qubit): gate_2 = getattr(gates, gate_2)(qubit) initial_state = random_statevector(2 ** (qubit + 1), backend=backend) final_state = apply_gates(backend, [gate_1, gate_2], qubit + 1, initial_state) - backend.assert_allclose(final_state, initial_state) + backend.assert_allclose(final_state, initial_state, atol=1e-6) @pytest.mark.parametrize("nqubits", [1, 2]) def test_unitary_dagger(backend, nqubits): matrix = np.random.random((2**nqubits, 2**nqubits)) + matrix = backend.cast(matrix) gate = gates.Unitary(matrix, *range(nqubits)) initial_state = random_statevector(2**nqubits, backend=backend) final_state = apply_gates(backend, [gate, gate.dagger()], nqubits, initial_state) - matrix = backend.cast(matrix, dtype=matrix.dtype) - target_state = np.dot(matrix, initial_state) - target_state = np.dot(np.conj(matrix).T, target_state) - backend.assert_allclose(final_state, target_state) + target_state = backend.np.matmul(matrix, initial_state) + target_state = backend.np.matmul(backend.np.conj(matrix).T, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_controlled_unitary_dagger(backend): @@ -1670,11 +1673,12 @@ def test_controlled_unitary_dagger(backend): matrix = np.random.random((2, 2)) matrix = expm(1j * (matrix + matrix.T)) + matrix = backend.cast(matrix) gate = gates.Unitary(matrix, 0).controlled_by(1, 2, 3, 4) nqubits = 5 initial_state = random_statevector(2**nqubits, backend=backend) final_state = apply_gates(backend, [gate, gate.dagger()], 5, initial_state) - backend.assert_allclose(final_state, initial_state) + backend.assert_allclose(final_state, initial_state, atol=1e-6) def test_generalizedfsim_dagger(backend): @@ -1687,7 +1691,7 @@ def test_generalizedfsim_dagger(backend): nqubits = 2 initial_state = random_statevector(2**nqubits, backend=backend) final_state = apply_gates(backend, [gate, gate.dagger()], 2, initial_state) - backend.assert_allclose(final_state, initial_state) + backend.assert_allclose(final_state, initial_state, atol=1e-6) ############################################################################### @@ -1728,9 +1732,9 @@ def test_x_decomposition_execution(backend, target, controls, free, use_toffolis gate = gates.X(target).controlled_by(*controls) nqubits = max((target,) + controls + free) + 1 initial_state = random_statevector(2**nqubits, backend=backend) - target_state = backend.apply_gate(gate, np.copy(initial_state), nqubits) + target_state = backend.apply_gate(gate, backend.np.copy(initial_state), nqubits) dgates = gate.decompose(*free, use_toffolis=use_toffolis) - final_state = np.copy(initial_state) + final_state = backend.np.copy(initial_state) for gate in dgates: final_state = backend.apply_gate(gate, final_state, nqubits) backend.assert_allclose(final_state, target_state, atol=1e-6) diff --git a/tests/test_hamiltonians.py b/tests/test_hamiltonians.py index c0089cc59b..5dce386185 100644 --- a/tests/test_hamiltonians.py +++ b/tests/test_hamiltonians.py @@ -456,5 +456,5 @@ def test_hamiltonian_energy_fluctuation(backend): gs_energy_fluctuation = ham.energy_fluctuation(ground_state) zs_energy_fluctuation = ham.energy_fluctuation(zero_state) - assert np.isclose(gs_energy_fluctuation, 0, atol=1e-5) + assert np.isclose(backend.to_numpy(gs_energy_fluctuation), 0, atol=1e-5) assert gs_energy_fluctuation < zs_energy_fluctuation diff --git a/tests/test_hamiltonians_terms.py b/tests/test_hamiltonians_terms.py index 1678e6d26e..8cf6f10771 100644 --- a/tests/test_hamiltonians_terms.py +++ b/tests/test_hamiltonians_terms.py @@ -48,10 +48,10 @@ def test_hamiltonian_term_gates(backend): backend.assert_allclose(gate.matrix(backend), matrix) initial_state = random_statevector(2**nqubits, backend=backend) - final_state = term(backend, np.copy(initial_state), nqubits) + final_state = term(backend, backend.np.copy(initial_state), nqubits) c = models.Circuit(nqubits) c.add(gates.Unitary(matrix, 1, 2)) - target_state = backend.execute_circuit(c, np.copy(initial_state)).state() + target_state = backend.execute_circuit(c, backend.np.copy(initial_state)).state() backend.assert_allclose(final_state, target_state) @@ -65,9 +65,9 @@ def test_hamiltonian_term_exponentiation(backend): backend.assert_allclose(term.exp(0.5), exp_matrix) initial_state = random_statevector(4, backend=backend) - final_state = term(backend, np.copy(initial_state), 2, term.expgate(0.5)) + final_state = term(backend, backend.np.copy(initial_state), 2, term.expgate(0.5)) exp_gate = gates.Unitary(exp_matrix, 1) - target_state = backend.apply_gate(exp_gate, np.copy(initial_state), 2) + target_state = backend.apply_gate(exp_gate, backend.np.copy(initial_state), 2) backend.assert_allclose(final_state, target_state) @@ -187,9 +187,9 @@ def test_symbolic_term_call(backend, density_matrix): else random_statevector(2**3, backend=backend) ) final_state = term( - backend, np.copy(initial_state), 3, density_matrix=density_matrix + backend, backend.np.copy(initial_state), 3, density_matrix=density_matrix ) - target_state = 2 * np.copy(initial_state) + target_state = 2 * backend.np.copy(initial_state) for matrix in matrixlist: target_state = matrix @ target_state backend.assert_allclose(final_state, target_state) diff --git a/tests/test_hamiltonians_trotter.py b/tests/test_hamiltonians_trotter.py index a69e9c4165..0c2d2f1756 100644 --- a/tests/test_hamiltonians_trotter.py +++ b/tests/test_hamiltonians_trotter.py @@ -123,23 +123,25 @@ def test_trotter_hamiltonian_three_qubit_term(backend): # Test that the `TrotterHamiltonian` dense matrix is correct eye = np.eye(2, dtype=complex) eye = backend.cast(eye, dtype=eye.dtype) - mm1 = np.kron(m1, eye) - mm2 = np.kron(np.kron(eye, eye), m2) - mm3 = np.kron(np.kron(eye, m3), np.kron(eye, eye)) + mm1 = backend.np.kron(m1, eye) + mm2 = backend.np.kron(backend.np.kron(eye, eye), m2) + mm3 = backend.np.kron(backend.np.kron(eye, m3), backend.np.kron(eye, eye)) target_ham = hamiltonians.Hamiltonian(4, mm1 + mm2 + mm3, backend=backend) backend.assert_allclose(ham.matrix, target_ham.matrix) dt = 1e-2 initial_state = random_statevector(2**4, backend=backend) circuit = ham.circuit(dt=dt) - final_state = backend.execute_circuit(circuit, np.copy(initial_state)).state() + final_state = backend.execute_circuit( + circuit, backend.np.copy(initial_state) + ).state() mm1 = backend.to_numpy(mm1) mm2 = backend.to_numpy(mm2) mm3 = backend.to_numpy(mm3) u = [expm(-0.5j * dt * (mm1 + mm3)), expm(-0.5j * dt * mm2)] u = backend.cast(u) - target_state = np.dot(u[1], np.dot(u[0], initial_state)) - target_state = np.dot(u[0], np.dot(u[1], target_state)) + target_state = backend.np.matmul(u[1], backend.np.matmul(u[0], initial_state)) + target_state = backend.np.matmul(u[0], backend.np.matmul(u[1], target_state)) backend.assert_allclose(final_state, target_state) @@ -150,7 +152,7 @@ def test_symbolic_hamiltonian_circuit_different_dts(backend): b = ham.circuit(0.1) matrix1 = ham.circuit(0.2).unitary(backend) matrix2 = (a + b).unitary(backend) - np.testing.assert_allclose(matrix1, matrix2) + backend.assert_allclose(matrix1, matrix2) def test_old_trotter_hamiltonian_errors(): diff --git a/tests/test_measurements_collapse.py b/tests/test_measurements_collapse.py index e0422bde8e..34cc48d9c7 100644 --- a/tests/test_measurements_collapse.py +++ b/tests/test_measurements_collapse.py @@ -41,10 +41,10 @@ def assign_value(rho, index, value): initial_rho = random_density_matrix(2**nqubits, backend=backend) c = models.Circuit(nqubits, density_matrix=True) r = c.add(gates.M(*targets, collapse=True)) - final_rho = backend.execute_circuit(c, np.copy(initial_rho), nshots=1) + final_rho = backend.execute_circuit(c, backend.np.copy(initial_rho), nshots=1) samples = r.samples()[0] - target_rho = np.reshape(initial_rho, 2 * nqubits * (2,)) + target_rho = backend.np.reshape(initial_rho, 2 * nqubits * (2,)) for q, r in zip(targets, samples): r = int(r) slicer = 2 * nqubits * [slice(None)] @@ -54,8 +54,8 @@ def assign_value(rho, index, value): target_rho = assign_value(target_rho, tuple(slicer), 0) slicer[q], slicer[q + nqubits] = 1 - r, r target_rho = assign_value(target_rho, tuple(slicer), 0) - target_rho = np.reshape(target_rho, initial_rho.shape) - target_rho = target_rho / np.trace(target_rho) + target_rho = backend.np.reshape(target_rho, initial_rho.shape) + target_rho = target_rho / backend.np.trace(target_rho) backend.assert_allclose(final_rho, target_rho) @@ -99,14 +99,14 @@ def test_measurement_result_parameters_random(backend): c.add(gates.RY(0, theta=np.pi * r.symbols[0] / 5)) c.add(gates.RX(2, theta=np.pi * r.symbols[0] / 4)) final_state = backend.execute_circuit( - c, initial_state=np.copy(initial_state), nshots=1 + c, initial_state=backend.np.copy(initial_state), nshots=1 ) backend.set_seed(123) c = models.Circuit(4, density_matrix=True) m = c.add(gates.M(1, collapse=True)) target_state = backend.execute_circuit( - c, initial_state=np.copy(initial_state), nshots=1 + c, initial_state=backend.np.copy(initial_state), nshots=1 ).state() if int(m.symbols[0].outcome()): c = models.Circuit(4, density_matrix=True) @@ -127,13 +127,13 @@ def test_measurement_result_parameters_repeated_execution(backend, use_loop): final_states = [] for _ in range(20): final_state = backend.execute_circuit( - c, initial_state=np.copy(initial_state), nshots=1 + c, initial_state=backend.np.copy(initial_state), nshots=1 ) final_states.append(final_state.state()) final_states = backend.np.mean(backend.cast(final_states), 0) else: final_states = backend.execute_circuit( - c, initial_state=np.copy(initial_state), nshots=20 + c, initial_state=backend.np.copy(initial_state), nshots=20 ).state() backend.set_seed(123) @@ -142,7 +142,7 @@ def test_measurement_result_parameters_repeated_execution(backend, use_loop): c = models.Circuit(4, density_matrix=True) m = c.add(gates.M(1, collapse=True)) target_state = backend.execute_circuit( - c, np.copy(initial_state), nshots=1 + c, backend.np.copy(initial_state), nshots=1 ).state() if int(m.symbols[0].outcome()): target_state = backend.apply_gate_density_matrix( @@ -173,7 +173,7 @@ def test_measurement_result_parameters_repeated_execution_final_measurements(bac c = models.Circuit(4, density_matrix=True) m = c.add(gates.M(1, collapse=True)) target_state = backend.execute_circuit( - c, np.copy(initial_state), nshots=1 + c, backend.np.copy(initial_state), nshots=1 ).state() c = models.Circuit(4, density_matrix=True) if int(m.symbols[0].outcome()): @@ -192,12 +192,14 @@ def test_measurement_result_parameters_multiple_qubits(backend): r = c.add(gates.M(0, 1, 2, collapse=True)) c.add(gates.RY(1, theta=np.pi * r.symbols[0] / 5)) c.add(gates.RX(3, theta=np.pi * r.symbols[2] / 3)) - final_state = backend.execute_circuit(c, np.copy(initial_state), nshots=1) + final_state = backend.execute_circuit(c, backend.np.copy(initial_state), nshots=1) backend.set_seed(123) c = models.Circuit(4, density_matrix=True) m = c.add(gates.M(0, 1, 2, collapse=True)) - target_state = backend.execute_circuit(c, np.copy(initial_state), nshots=1).state() + target_state = backend.execute_circuit( + c, backend.np.copy(initial_state), nshots=1 + ).state() # not including in coverage because outcomes are probabilistic and may # not occur for the CI run if int(m.symbols[0].outcome()): # pragma: no cover diff --git a/tests/test_models_circuit_execution.py b/tests/test_models_circuit_execution.py index 0158cf1001..297b715a1e 100644 --- a/tests/test_models_circuit_execution.py +++ b/tests/test_models_circuit_execution.py @@ -99,7 +99,7 @@ def test_density_matrix_circuit(backend): c.add(gates.H(1)) c.add(gates.CNOT(0, 1)) c.add(gates.H(2)) - final_rho = backend.execute_circuit(c, np.copy(initial_rho)).state() + final_rho = backend.execute_circuit(c, backend.np.copy(initial_rho)).state() h = np.array([[1, 1], [1, -1]]) / np.sqrt(2) cnot = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) @@ -107,11 +107,9 @@ def test_density_matrix_circuit(backend): m2 = np.kron(cnot, np.eye(2)) m3 = np.kron(np.eye(4), h) - m1 = backend.cast(m1, dtype=m1.dtype) - m2 = backend.cast(m2, dtype=m2.dtype) - m3 = backend.cast(m3, dtype=m3.dtype) - - target_rho = np.dot(m1, np.dot(initial_rho, np.transpose(np.conj(m1)))) + target_rho = np.dot( + m1, np.dot(backend.to_numpy(initial_rho), np.transpose(np.conj(m1))) + ) target_rho = np.dot(m2, np.dot(target_rho, np.transpose(np.conj(m2)))) target_rho = np.dot(m3, np.dot(target_rho, np.transpose(np.conj(m3)))) diff --git a/tests/test_models_circuit_features.py b/tests/test_models_circuit_features.py index b8a1ade5e7..c0b08f8c5d 100644 --- a/tests/test_models_circuit_features.py +++ b/tests/test_models_circuit_features.py @@ -159,10 +159,10 @@ def test_inverse_circuit_execution(backend, accelerators, fuse): else: c = c.fuse() invc = c.invert() - target_state = np.ones(2**4) / 4 + target_state = np.ones(2**4) / 4.0 final_state = backend.execute_circuit(c, initial_state=np.copy(target_state))._state final_state = backend.execute_circuit(invc, initial_state=final_state)._state - backend.assert_allclose(final_state, target_state) + backend.assert_allclose(final_state, target_state, atol=1e-6) def test_circuit_invert_and_addition_execution(backend, accelerators): diff --git a/tests/test_models_circuit_parametrized.py b/tests/test_models_circuit_parametrized.py index 6bc1698321..1a761f16eb 100644 --- a/tests/test_models_circuit_parametrized.py +++ b/tests/test_models_circuit_parametrized.py @@ -5,7 +5,6 @@ import numpy as np import pytest -import qibo from qibo import Circuit, gates @@ -21,14 +20,14 @@ def exact_state(theta): theta = 0.1234 gate = gates.RX(0, theta=theta) - initial_state = np.ones(2) / np.sqrt(2) + initial_state = backend.cast(np.ones(2) / np.sqrt(2)) final_state = backend.apply_gate(gate, initial_state, 1) target_state = exact_state(theta) backend.assert_allclose(final_state, target_state) theta = 0.4321 gate.parameters = theta - initial_state = np.ones(2) / np.sqrt(2) + initial_state = backend.cast(np.ones(2) / np.sqrt(2)) final_state = backend.apply_gate(gate, initial_state, 1) target_state = exact_state(theta) backend.assert_allclose(final_state, target_state) diff --git a/tests/test_models_dbi.py b/tests/test_models_dbi.py index 3b5d32664b..d19168caa9 100644 --- a/tests/test_models_dbi.py +++ b/tests/test_models_dbi.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from qibo import hamiltonians, set_backend +from qibo import hamiltonians from qibo.hamiltonians import Hamiltonian from qibo.models.dbi.double_bracket import ( DoubleBracketCostFunction, @@ -229,12 +229,14 @@ def test_params_to_diagonal_operator(backend, step): nqubits=nqubits, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict, + backend=backend, ), ) operator_element = params_to_diagonal_operator( params, nqubits=nqubits, parameterization=ParameterizationTypes.computational, + backend=backend, ) for i in range(len(params)): backend.assert_allclose( @@ -261,7 +263,7 @@ def test_gradient_descent(backend, order): pauli_operators = list(pauli_operator_dict.values()) # let initial d be approximation of $\Delta(H) d_coef_pauli = decompose_into_pauli_basis( - dbi.diagonal_h_matrix, pauli_operators=pauli_operators + dbi.diagonal_h_matrix, pauli_operators=pauli_operators, backend=backend ) d_pauli = sum([d_coef_pauli[i] * pauli_operators[i] for i in range(nqubits)]) loss_hist_pauli, d_params_hist_pauli, s_hist_pauli = gradient_descent( @@ -271,6 +273,7 @@ def test_gradient_descent(backend, order): ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict, pauli_parameterization_order=order, + backend=backend, ) assert loss_hist_pauli[-1] < initial_off_diagonal_norm @@ -281,6 +284,10 @@ def test_gradient_descent(backend, order): _, _, ) = gradient_descent( - dbi, NSTEPS, d_coef_computational_partial, ParameterizationTypes.computational + dbi, + NSTEPS, + d_coef_computational_partial, + ParameterizationTypes.computational, + backend=backend, ) assert loss_hist_computational_partial[-1] < initial_off_diagonal_norm diff --git a/tests/test_models_dbi_utils.py b/tests/test_models_dbi_utils.py new file mode 100644 index 0000000000..29d1b1058b --- /dev/null +++ b/tests/test_models_dbi_utils.py @@ -0,0 +1,61 @@ +""""Testing utils for DoubleBracketIteration model""" + +import numpy as np +import pytest + +from qibo.hamiltonians import Hamiltonian +from qibo.models.dbi.double_bracket import ( + DoubleBracketGeneratorType, + DoubleBracketIteration, +) +from qibo.models.dbi.utils import * +from qibo.models.dbi.utils_dbr_strategies import select_best_dbr_generator +from qibo.quantum_info import random_hermitian + +NSTEPS = 5 +"""Number of steps for evolution.""" + + +@pytest.mark.parametrize("nqubits", [1, 2]) +def test_generate_Z_operators(backend, nqubits): + h0 = random_hermitian(2**nqubits, backend=backend) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits=nqubits, matrix=h0, backend=backend) + ) + generate_Z = generate_Z_operators(nqubits, backend=backend) + Z_ops = list(generate_Z.values()) + + delta_h0 = dbi.diagonal_h_matrix + dephasing_channel = (sum([Z_op @ h0 @ Z_op for Z_op in Z_ops]) + h0) / 2**nqubits + norm_diff = np.linalg.norm(backend.to_numpy(delta_h0 - dephasing_channel)) + + assert norm_diff < 1e-3 + + +@pytest.mark.parametrize("nqubits", [1, 2]) +@pytest.mark.parametrize("step", [0.1, 0.2]) +def test_select_best_dbr_generator(backend, nqubits, step): + h0 = random_hermitian(2**nqubits, seed=1, backend=backend) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + ) + generate_Z = generate_Z_operators(nqubits, backend=backend) + Z_ops = list(generate_Z.values()) + initial_off_diagonal_norm = dbi.off_diagonal_norm + + for _ in range(NSTEPS): + dbi, idx, step_optimize, flip = select_best_dbr_generator( + dbi, Z_ops, step=step, compare_canonical=True, max_evals=5 + ) + + assert initial_off_diagonal_norm > dbi.off_diagonal_norm + + +def test_copy_dbi(backend): + h0 = random_hermitian(4, seed=1, backend=backend) + dbi = DoubleBracketIteration(Hamiltonian(2, h0, backend=backend)) + dbi_copy = copy_dbi_object(dbi) + + assert dbi is not dbi_copy + assert dbi.h.nqubits == dbi_copy.h.nqubits diff --git a/tests/test_models_encodings.py b/tests/test_models_encodings.py index cd8994df3e..1a935d75e3 100644 --- a/tests/test_models_encodings.py +++ b/tests/test_models_encodings.py @@ -132,7 +132,7 @@ def test_unary_encoder(backend, nqubits, architecture, kind): circuit = unary_encoder(data, architecture=architecture) state = backend.execute_circuit(circuit).state() - indexes = np.flatnonzero(state) + indexes = np.flatnonzero(backend.to_numpy(state)) state = backend.np.real(state[indexes]) backend.assert_allclose(state, backend.cast(data) / backend.calculate_norm(data, 2)) @@ -165,7 +165,7 @@ def test_unary_encoder_random_gaussian(backend, nqubits, seed): for _ in range(samples): circuit = unary_encoder_random_gaussian(nqubits, seed=local_state) state = backend.execute_circuit(circuit).state() - indexes = np.flatnonzero(state) + indexes = np.flatnonzero(backend.to_numpy(state)) state = np.real(state[indexes]) amplitudes += [float(elem) for elem in list(state)] diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py index 3881b2b0db..4335f909f7 100644 --- a/tests/test_models_error_mitigation.py +++ b/tests/test_models_error_mitigation.py @@ -133,7 +133,7 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout): backend=backend, ) - assert np.abs(exact - estimate) <= np.abs(exact - noisy) + assert backend.np.abs(exact - estimate) <= backend.np.abs(exact - noisy) @pytest.mark.parametrize("nqubits", [3]) @@ -192,7 +192,7 @@ def test_cdr(backend, nqubits, noise, full_output, readout): if full_output: estimate = estimate[0] - assert np.abs(exact - estimate) <= np.abs(exact - noisy) + assert backend.np.abs(exact - estimate) <= backend.np.abs(exact - noisy) @pytest.mark.parametrize("nqubits", [3]) @@ -275,7 +275,7 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout): if full_output: estimate = estimate[0] - assert np.abs(exact - estimate) <= np.abs(exact - noisy) + assert backend.np.abs(exact - estimate) <= backend.np.abs(exact - noisy) @pytest.mark.parametrize("nqubits,nmeas", [(3, 2)]) @@ -372,4 +372,4 @@ def test_ics(backend, nqubits, noise, full_output, readout): if full_output: estimate = estimate[0] - assert np.abs(exact - estimate) <= np.abs(exact - noisy) + assert backend.np.abs(exact - estimate) <= backend.np.abs(exact - noisy) diff --git a/tests/test_models_qft.py b/tests/test_models_qft.py index 0165b1a933..33a192adf6 100644 --- a/tests/test_models_qft.py +++ b/tests/test_models_qft.py @@ -55,7 +55,7 @@ def test_qft_execution(backend, accelerators, nqubits, random): initial_state = random_statevector(2**nqubits, backend=backend) else: initial_state = backend.zero_state(nqubits) - final_state = backend.execute_circuit(c, np.copy(initial_state))._state + final_state = backend.execute_circuit(c, backend.np.copy(initial_state))._state target_state = exact_qft(backend.to_numpy(initial_state)) backend.assert_allclose(final_state, target_state) diff --git a/tests/test_models_variational.py b/tests/test_models_variational.py index 2ee1d2bc59..c596990488 100644 --- a/tests/test_models_variational.py +++ b/tests/test_models_variational.py @@ -101,12 +101,10 @@ def myloss(parameters, circuit, target): @pytest.mark.parametrize(test_names, test_values) def test_vqe(backend, method, options, compile, filename): """Performs a VQE circuit minimization test.""" - if backend.name == "pytorch": - pytest.skip("Skipping VQE test for pytorch backend.") - if (method == "sgd" or compile) and backend.name != "tensorflow": + if (method == "sgd" or compile) and (backend.name not in ("tensorflow", "pytorch")): pytest.skip("Skipping SGD test for unsupported backend.") - if method != "sgd" and backend.name == "tensorflow": - pytest.skip("Skipping scipy optimizers for tensorflow.") + if method != "sgd" and backend.name in ("tensorflow", "pytorch"): + pytest.skip("Skipping scipy optimizers for tensorflow and pytorch.") n_threads = backend.nthreads backend.set_threads(1) nqubits = 3 @@ -127,6 +125,9 @@ def test_vqe(backend, method, options, compile, filename): hamiltonian = hamiltonians.XXZ(nqubits=nqubits, backend=backend) np.random.seed(0) initial_parameters = np.random.uniform(0, 2 * np.pi, 2 * nqubits * layers + nqubits) + initial_parameters = backend.cast( + initial_parameters, dtype=initial_parameters.dtype + ) v = models.VQE(circuit, hamiltonian) loss_values = [] @@ -183,7 +184,7 @@ def test_qaoa_execution(backend, solver, dense, accel=None): else: atol = 0 - target_state = np.copy(state) + target_state = backend.cast(state, copy=True) h_matrix = backend.to_numpy(h.matrix) m_matrix = backend.to_numpy(m.matrix) for i, p in enumerate(params): @@ -268,11 +269,14 @@ def test_qaoa_errors(backend): @pytest.mark.parametrize(test_names, test_values) def test_qaoa_optimization(backend, method, options, dense, filename): - if method == "sgd" and backend.name != "tensorflow": + if (method == "sgd") and (backend.name not in ["tensorflow", "pytorch"]): pytest.skip("Skipping SGD test for unsupported backend.") + if method != "sgd" and backend.name in ("tensorflow", "pytorch"): + pytest.skip("Skipping scipy optimizers for tensorflow and pytorch.") h = hamiltonians.XXZ(3, dense=dense, backend=backend) qaoa = models.QAOA(h) initial_p = [0.05, 0.06, 0.07, 0.08] + initial_p = backend.cast(initial_p, dtype=np.float64) best, params, _ = qaoa.minimize(initial_p, method=method, options=options) if filename is not None: assert_regression_fixture(backend, params, filename) @@ -319,8 +323,7 @@ def __call__(self, x): @pytest.mark.parametrize(test_names, test_values) def test_aavqe(backend, method, options, compile, filename): """Performs a AAVQE circuit minimization test.""" - if backend.name == "pytorch": - pytest.skip("Skipping VQE test for pytorch backend.") + nqubits = 4 layers = 1 circuit = models.Circuit(nqubits) diff --git a/tests/test_noise.py b/tests/test_noise.py index a7882050d3..325ec95e2d 100644 --- a/tests/test_noise.py +++ b/tests/test_noise.py @@ -67,12 +67,12 @@ def test_kraus_error(backend, density_matrix, nshots): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), initial_state=np.copy(initial_psi), nshots=nshots + noise.apply(circuit), initial_state=backend.np.copy(initial_psi), nshots=nshots ) final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -125,12 +125,12 @@ def test_unitary_error(backend, density_matrix, nshots): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), initial_state=np.copy(initial_psi), nshots=nshots + noise.apply(circuit), initial_state=backend.np.copy(initial_psi), nshots=nshots ) final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -180,12 +180,12 @@ def test_pauli_error(backend, density_matrix, nshots): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), initial_state=np.copy(initial_psi), nshots=nshots + noise.apply(circuit), initial_state=backend.np.copy(initial_psi), nshots=nshots ) final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -237,7 +237,7 @@ def test_depolarizing_error(backend, density_matrix, nshots): final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -282,11 +282,11 @@ def test_thermal_error(backend, density_matrix): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), np.copy(initial_psi) + noise.apply(circuit), backend.np.copy(initial_psi) )._state backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, np.copy(initial_psi) + target_circuit, backend.np.copy(initial_psi) )._state backend.assert_allclose(final_state, target_final_state) @@ -335,7 +335,7 @@ def test_amplitude_damping_error(backend, density_matrix, nshots): final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -388,7 +388,7 @@ def test_phase_damping_error(backend, density_matrix, nshots): final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -416,11 +416,11 @@ def test_readout_error(backend, density_matrix): circuit = Circuit(nqubits, density_matrix=density_matrix) circuit.add(gates.M(0)) final_state = backend.execute_circuit( - noise.apply(circuit), initial_state=np.copy(state) + noise.apply(circuit), initial_state=backend.np.copy(state) ) target_state = gates.ReadoutErrorChannel(0, P).apply_density_matrix( - backend, np.copy(state), nqubits + backend, backend.np.copy(state), nqubits ) backend.assert_allclose(final_state, target_state) @@ -454,11 +454,11 @@ def test_reset_error(backend, density_matrix): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), np.copy(initial_psi) + noise.apply(circuit), backend.np.copy(initial_psi) )._state backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, np.copy(initial_psi) + target_circuit, backend.np.copy(initial_psi) )._state backend.assert_allclose(final_state, target_final_state) @@ -505,12 +505,12 @@ def test_custom_error(backend, density_matrix, nshots): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), initial_state=np.copy(initial_psi), nshots=nshots + noise.apply(circuit), initial_state=backend.np.copy(initial_psi), nshots=nshots ) final_state_samples = final_state.samples() if nshots else None backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi), nshots=nshots + target_circuit, initial_state=backend.np.copy(initial_psi), nshots=nshots ) target_final_state_samples = target_final_state.samples() if nshots else None @@ -564,11 +564,11 @@ def condition_3_pi_2(gate): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), np.copy(initial_psi) + noise.apply(circuit), backend.np.copy(initial_psi) )._state backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, np.copy(initial_psi) + target_circuit, backend.np.copy(initial_psi) )._state backend.assert_allclose(final_state, target_final_state) @@ -606,11 +606,11 @@ def test_gate_independent_noise(backend, density_matrix): ) backend.set_seed(123) final_state = backend.execute_circuit( - noise.apply(circuit), initial_state=np.copy(initial_psi) + noise.apply(circuit), initial_state=backend.np.copy(initial_psi) )._state backend.set_seed(123) target_final_state = backend.execute_circuit( - target_circuit, initial_state=np.copy(initial_psi) + target_circuit, initial_state=backend.np.copy(initial_psi) )._state backend.assert_allclose(final_state, target_final_state) diff --git a/tests/test_quantum_info_basis.py b/tests/test_quantum_info_basis.py index bbca32bc68..6d35863504 100644 --- a/tests/test_quantum_info_basis.py +++ b/tests/test_quantum_info_basis.py @@ -53,18 +53,19 @@ def test_pauli_basis( if vectorize: basis_test = [ - vectorization(matrix, order=order, backend=backend) for matrix in basis_test + vectorization(backend.cast(matrix), order=order, backend=backend) + for matrix in basis_test ] basis_test = backend.cast(basis_test) if normalize: - basis_test /= np.sqrt(2**nqubits) + basis_test = basis_test / np.sqrt(2**nqubits) if vectorize and sparse: elements, indexes = [], [] for row in basis_test: - row_indexes = list(np.flatnonzero(row)) + row_indexes = list(np.flatnonzero(backend.to_numpy(row))) indexes.append(row_indexes) elements.append(row[row_indexes]) indexes = backend.cast(indexes) @@ -82,15 +83,19 @@ def test_pauli_basis( elements, indexes, basis[0], basis[1] ): backend.assert_allclose( - np.linalg.norm(elem_test - elem) < PRECISION_TOL, True + np.linalg.norm(backend.to_numpy(elem_test - elem)) < PRECISION_TOL, + True, ) backend.assert_allclose( - np.linalg.norm(ind_test - ind) < PRECISION_TOL, True + np.linalg.norm(backend.to_numpy(ind_test - ind)) < PRECISION_TOL, + True, ) else: for pauli, pauli_test in zip(basis, basis_test): backend.assert_allclose( - np.linalg.norm(pauli - pauli_test) < PRECISION_TOL, True + np.linalg.norm(backend.to_numpy(pauli - pauli_test)) + < PRECISION_TOL, + True, ) comp_basis_to_pauli(nqubits, normalize, sparse, order, pauli_order, backend) diff --git a/tests/test_quantum_info_entanglement.py b/tests/test_quantum_info_entanglement.py index a98cf1909b..a641d01c94 100644 --- a/tests/test_quantum_info_entanglement.py +++ b/tests/test_quantum_info_entanglement.py @@ -50,7 +50,7 @@ def test_concurrence_and_formation(backend, bipartition, base, check_purity): backend.assert_allclose(0.0 <= concur <= np.sqrt(2), True) backend.assert_allclose(0.0 <= ent_form <= 1.0, True) - state = np.kron( + state = backend.np.kron( random_density_matrix(2, pure=True, backend=backend), random_density_matrix(2, pure=True, backend=backend), ) diff --git a/tests/test_quantum_info_entropies.py b/tests/test_quantum_info_entropies.py index 5013f4f70d..84129e3c6d 100644 --- a/tests/test_quantum_info_entropies.py +++ b/tests/test_quantum_info_entropies.py @@ -524,9 +524,11 @@ def test_renyi_entropy(backend, alpha, base): target = von_neumann_entropy(state, base=base, backend=backend) elif alpha == np.inf: target = backend.calculate_norm_density_matrix(state, order=2) - target = -1 * np.log2(target) / np.log2(base) + target = -1 * backend.np.log2(target) / np.log2(base) else: - target = np.log2(np.trace(np.linalg.matrix_power(state, alpha))) + target = np.log2( + np.trace(np.linalg.matrix_power(backend.to_numpy(state), alpha)) + ) target = (1 / (1 - alpha)) * target / np.log2(base) backend.assert_allclose( @@ -605,7 +607,7 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag): new_state = _matrix_power(state, 0.5, backend) new_target = _matrix_power(target, 0.5, backend) - log = np.log2( + log = backend.np.log2( backend.calculate_norm_density_matrix( new_state @ new_target, order=1 ) @@ -614,14 +616,14 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag): log = -2 * log / np.log2(base) else: if len(state.shape) == 1: - state = np.outer(state, np.conj(state)) + state = backend.np.outer(state, backend.np.conj(state)) if len(target.shape) == 1: - target = np.outer(target, np.conj(target)) + target = backend.np.outer(target, backend.np.conj(target)) log = _matrix_power(state, alpha, backend) log = log @ _matrix_power(target, 1 - alpha, backend) - log = np.log2(np.trace(log)) + log = backend.np.log2(backend.np.trace(log)) log = (1 / (alpha - 1)) * log / np.log2(base) @@ -666,7 +668,7 @@ def test_tsallis_entropy(backend, alpha, base): target = von_neumann_entropy(state, base=base, backend=backend) else: target = (1 / (1 - alpha)) * ( - np.trace(_matrix_power(state, alpha, backend)) - 1 + backend.np.trace(_matrix_power(state, alpha, backend)) - 1 ) backend.assert_allclose( @@ -741,7 +743,7 @@ def test_entanglement_entropy(backend, bipartition, base, check_hermitian): backend.assert_allclose(entang_entrop, test, atol=PRECISION_TOL) # Product state - state = np.kron( + state = backend.np.kron( random_statevector(2, backend=backend), random_statevector(2, backend=backend) ) @@ -752,5 +754,4 @@ def test_entanglement_entropy(backend, bipartition, base, check_hermitian): check_hermitian=check_hermitian, backend=backend, ) - backend.assert_allclose(entang_entrop, 0.0, atol=PRECISION_TOL) diff --git a/tests/test_quantum_info_metrics.py b/tests/test_quantum_info_metrics.py index ee0dda211f..b5b28c60bf 100644 --- a/tests/test_quantum_info_metrics.py +++ b/tests/test_quantum_info_metrics.py @@ -32,23 +32,27 @@ def test_purity_and_impurity(backend): with pytest.raises(TypeError): state = np.random.rand(2, 3) state = backend.cast(state, dtype=state.dtype) - test = purity(state) + test = purity(state, backend=backend) state = np.array([1.0, 0.0, 0.0, 0.0]) state = backend.cast(state, dtype=state.dtype) - backend.assert_allclose(purity(state), 1.0, atol=PRECISION_TOL) - backend.assert_allclose(impurity(state), 0.0, atol=PRECISION_TOL) + backend.assert_allclose(purity(state, backend=backend), 1.0, atol=PRECISION_TOL) + backend.assert_allclose(impurity(state, backend=backend), 0.0, atol=PRECISION_TOL) - state = np.outer(np.conj(state), state) + state = backend.np.outer(backend.np.conj(state), state) state = backend.cast(state, dtype=state.dtype) - backend.assert_allclose(purity(state), 1.0, atol=PRECISION_TOL) - backend.assert_allclose(impurity(state), 0.0, atol=PRECISION_TOL) + backend.assert_allclose(purity(state, backend=backend), 1.0, atol=PRECISION_TOL) + backend.assert_allclose(impurity(state, backend=backend), 0.0, atol=PRECISION_TOL) dim = 4 state = backend.identity_density_matrix(2) state = backend.cast(state, dtype=state.dtype) - backend.assert_allclose(purity(state), 1.0 / dim, atol=PRECISION_TOL) - backend.assert_allclose(impurity(state), 1.0 - 1.0 / dim, atol=PRECISION_TOL) + backend.assert_allclose( + purity(state, backend=backend), 1.0 / dim, atol=PRECISION_TOL + ) + backend.assert_allclose( + impurity(state, backend=backend), 1.0 - 1.0 / dim, atol=PRECISION_TOL + ) @pytest.mark.parametrize("check_hermitian", [False, True]) @@ -90,10 +94,8 @@ def test_trace_distance(backend, check_hermitian): atol=PRECISION_TOL, ) - state = np.outer(np.conj(state), state) - target = np.outer(np.conj(target), target) - state = backend.cast(state, dtype=state.dtype) - target = backend.cast(target, dtype=target.dtype) + state = backend.np.outer(backend.np.conj(state), state) + target = backend.np.outer(backend.np.conj(target), target) backend.assert_allclose( trace_distance(state, target, check_hermitian=check_hermitian, backend=backend), 0.0, @@ -136,19 +138,23 @@ def test_hilbert_schmidt_distance(backend): target = np.array([1.0, 0.0, 0.0, 0.0]) state = backend.cast(state, dtype=state.dtype) target = backend.cast(target, dtype=target.dtype) - backend.assert_allclose(hilbert_schmidt_distance(state, target), 0.0) + backend.assert_allclose( + hilbert_schmidt_distance(state, target, backend=backend), 0.0 + ) - state = np.outer(np.conj(state), state) - target = np.outer(np.conj(target), target) - state = backend.cast(state, dtype=state.dtype) - target = backend.cast(target, dtype=target.dtype) - backend.assert_allclose(hilbert_schmidt_distance(state, target), 0.0) + state = backend.np.outer(backend.np.conj(state), state) + target = backend.np.outer(backend.np.conj(target), target) + backend.assert_allclose( + hilbert_schmidt_distance(state, target, backend=backend), 0.0 + ) state = np.array([0.0, 1.0, 0.0, 0.0]) target = np.array([1.0, 0.0, 0.0, 0.0]) state = backend.cast(state, dtype=state.dtype) target = backend.cast(target, dtype=target.dtype) - backend.assert_allclose(hilbert_schmidt_distance(state, target), 2.0) + backend.assert_allclose( + hilbert_schmidt_distance(state, target, backend=backend), 2.0 + ) @pytest.mark.parametrize("check_hermitian", [True, False]) @@ -201,10 +207,8 @@ def test_fidelity_and_infidelity_and_bures(backend, check_hermitian): atol=PRECISION_TOL, ) - state = np.outer(np.conj(state), state) - target = np.outer(np.conj(target), target) - state = backend.cast(state, dtype=state.dtype) - target = backend.cast(target, dtype=target.dtype) + state = backend.np.outer(backend.np.conj(state), state) + target = backend.np.outer(backend.np.conj(target), target) backend.assert_allclose( fidelity(state, target, check_hermitian, backend=backend), 1.0, diff --git a/tests/test_quantum_info_quantum_networks.py b/tests/test_quantum_info_quantum_networks.py index 1f20a97b45..5ebbbea3d3 100644 --- a/tests/test_quantum_info_quantum_networks.py +++ b/tests/test_quantum_info_quantum_networks.py @@ -273,7 +273,7 @@ def test_apply(backend): network = QuantumNetwork(unitary, (dims, dims), pure=True, backend=backend) applied = network.apply(state) - target = unitary @ state @ np.transpose(np.conj(unitary)) + target = unitary @ state @ backend.np.conj(unitary).T backend.assert_allclose(applied, target, atol=1e-8) diff --git a/tests/test_quantum_info_random.py b/tests/test_quantum_info_random.py index f974f57dd3..d522f1458d 100644 --- a/tests/test_quantum_info_random.py +++ b/tests/test_quantum_info_random.py @@ -49,9 +49,9 @@ def test_uniform_sampling_U3(backend, seed): expectation_values.append( [ - np.conj(state) @ X @ state, - np.conj(state) @ Y @ state, - np.conj(state) @ Z @ state, + backend.np.conj(state) @ X @ state, + backend.np.conj(state) @ Y @ state, + backend.np.conj(state) @ Z @ state, ] ) expectation_values = backend.cast(expectation_values) @@ -111,40 +111,40 @@ def test_random_hermitian(backend): # test if function returns Hermitian operator dims = 4 matrix = random_hermitian(dims, backend=backend) - matrix_dagger = np.transpose(np.conj(matrix)) + matrix_dagger = backend.np.conj(matrix).T norm = float(backend.calculate_norm_density_matrix(matrix - matrix_dagger, order=2)) backend.assert_allclose(norm < PRECISION_TOL, True) # test if function returns semidefinite Hermitian operator dims = 4 matrix = random_hermitian(dims, semidefinite=True, backend=backend) - matrix_dagger = np.transpose(np.conj(matrix)) + matrix_dagger = backend.np.conj(matrix).T norm = float(backend.calculate_norm_density_matrix(matrix - matrix_dagger, order=2)) backend.assert_allclose(norm < PRECISION_TOL, True) - eigenvalues = np.linalg.eigvalsh(matrix) + eigenvalues = np.linalg.eigvalsh(backend.to_numpy(matrix)) eigenvalues = np.real(eigenvalues) backend.assert_allclose(all(eigenvalues >= 0), True) # test if function returns normalized Hermitian operator dims = 4 matrix = random_hermitian(dims, normalize=True, backend=backend) - matrix_dagger = np.transpose(np.conj(matrix)) + matrix_dagger = backend.np.conj(matrix).T norm = float(backend.calculate_norm_density_matrix(matrix - matrix_dagger, order=2)) backend.assert_allclose(norm < PRECISION_TOL, True) - eigenvalues = np.linalg.eigvalsh(matrix) + eigenvalues = np.linalg.eigvalsh(backend.to_numpy(matrix)) eigenvalues = np.real(eigenvalues) backend.assert_allclose(all(eigenvalues <= 1), True) # test if function returns normalized and semidefinite Hermitian operator dims = 4 matrix = random_hermitian(dims, semidefinite=True, normalize=True, backend=backend) - matrix_dagger = np.transpose(np.conj(matrix)) + matrix_dagger = backend.np.conj(matrix).T norm = float(backend.calculate_norm(matrix - matrix_dagger, order=2)) backend.assert_allclose(norm < PRECISION_TOL, True) - eigenvalues = np.linalg.eigvalsh(matrix) + eigenvalues = np.linalg.eigvalsh(backend.to_numpy(matrix)) eigenvalues = np.real(eigenvalues) backend.assert_allclose(all(eigenvalues >= 0), True) backend.assert_allclose(all(eigenvalues <= 1), True) @@ -171,7 +171,7 @@ def test_random_unitary(backend, measure): # tests if operator is unitary (measure == "haar") dims = 4 matrix = random_unitary(dims, measure=measure, backend=backend) - matrix_dagger = np.transpose(np.conj(matrix)) + matrix_dagger = backend.np.conj(matrix).T matrix_inv = ( backend.np.inverse(matrix) if backend.name == "pytorch" @@ -234,7 +234,9 @@ def test_random_statevector(backend, seed): # tests if random statevector is a pure state dims = 4 state = random_statevector(dims, seed=seed, backend=backend) - backend.assert_allclose(abs(purity(state) - 1.0) < PRECISION_TOL, True) + backend.assert_allclose( + abs(purity(state, backend=backend) - 1.0) < PRECISION_TOL, True + ) @pytest.mark.parametrize("normalize", [False, True]) @@ -287,22 +289,31 @@ def test_random_density_matrix(backend, dims, pure, metric, basis, normalize): ) if basis is None and normalize is False: backend.assert_allclose( - np.real(np.trace(state)) <= 1.0 + PRECISION_TOL, True + np.real(np.trace(backend.to_numpy(state))) <= 1.0 + PRECISION_TOL, True + ) + backend.assert_allclose( + np.real(np.trace(backend.to_numpy(state))) >= 1.0 - PRECISION_TOL, True ) backend.assert_allclose( - np.real(np.trace(state)) >= 1.0 - PRECISION_TOL, True + purity(state, backend=backend) <= 1.0 + PRECISION_TOL, True ) - backend.assert_allclose(purity(state) <= 1.0 + PRECISION_TOL, True) if pure is True: - backend.assert_allclose(purity(state) >= 1.0 - PRECISION_TOL, True) - - state_dagger = np.transpose(np.conj(state)) - norm = float(norm_function(state - state_dagger, order=2)) + backend.assert_allclose( + purity(state, backend=backend) >= 1.0 - PRECISION_TOL, True + ) + norm = np.abs( + backend.to_numpy( + norm_function(state - backend.np.conj(state).T, order=2) + ) + ) backend.assert_allclose(norm < PRECISION_TOL, True) else: normalization = 1.0 if normalize is False else 1.0 / np.sqrt(dims) backend.assert_allclose(state[0], normalization) - assert all(np.abs(exp_value) <= normalization for exp_value in state[1:]) + assert all( + np.abs(backend.to_numpy(exp_value)) <= normalization + for exp_value in state[1:] + ) @pytest.mark.parametrize("seed", [10]) @@ -403,7 +414,11 @@ def test_pauli_single(backend): matrix = backend.cast(matrix, dtype=matrix.dtype) backend.assert_allclose( - float(backend.calculate_norm_density_matrix(matrix - result, order=2)) + np.abs( + backend.to_numpy( + backend.calculate_norm_density_matrix(matrix - result, order=2) + ) + ) < PRECISION_TOL, True, ) @@ -460,10 +475,9 @@ def test_random_pauli( True, ) else: - matrix = np.transpose(matrix, (1, 0, 2, 3)) + matrix = backend.np.transpose(matrix, (1, 0, 2, 3)) matrix = [reduce(backend.np.kron, row) for row in matrix] - dot = backend.np.matmul if backend.name == "pytorch" else np.dot - matrix = reduce(dot, matrix) + matrix = reduce(backend.np.matmul, matrix) if subset is None: backend.assert_allclose( @@ -516,10 +530,15 @@ def test_random_pauli_hamiltonian( ) if normalize is True: - backend.assert_allclose(np.abs(eigenvalues[0]) < PRECISION_TOL, True) - backend.assert_allclose(np.abs(eigenvalues[1] - 1) < PRECISION_TOL, True) backend.assert_allclose( - np.abs(eigenvalues[-1] - max_eigenvalue) < PRECISION_TOL, True + np.abs(backend.to_numpy(eigenvalues[0])) < PRECISION_TOL, True + ) + backend.assert_allclose( + np.abs(backend.to_numpy(eigenvalues[1]) - 1) < PRECISION_TOL, True + ) + backend.assert_allclose( + np.abs(backend.to_numpy(eigenvalues[-1]) - max_eigenvalue) < PRECISION_TOL, + True, ) diff --git a/tests/test_quantum_info_superoperator_transformations.py b/tests/test_quantum_info_superoperator_transformations.py index 38511a8a69..d93db2e5e6 100644 --- a/tests/test_quantum_info_superoperator_transformations.py +++ b/tests/test_quantum_info_superoperator_transformations.py @@ -47,14 +47,12 @@ @pytest.mark.parametrize("order", ["row", "column", "system"]) @pytest.mark.parametrize("nqubits", [1, 2, 3]) -def test_vectorization(backend, nqubits, order): +@pytest.mark.parametrize("statevector", [True, False]) +def test_vectorization(backend, nqubits, order, statevector): with pytest.raises(TypeError): - vectorization(np.array([[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]), backend=backend) - with pytest.raises(TypeError): - vectorization( - np.array([[[0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0]]], dtype="object"), - backend=backend, - ) + x = np.array([[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]) + x = backend.cast(x) + vectorization(x, backend=backend) with pytest.raises(TypeError): vectorization(np.array([]), backend=backend) with pytest.raises(TypeError): @@ -67,10 +65,16 @@ def test_vectorization(backend, nqubits, order): dim = 2**nqubits if nqubits == 1: - if order == "system" or order == "column": - matrix_test = [0, 2, 1, 3] + if statevector: + if order == "system": + matrix_test = [0, 0, 0, 1, 0, 0, 2, 3, 0, 2, 0, 3, 4, 6, 6, 9] + else: + matrix_test = [0, 0, 0, 0, 0, 1, 2, 3, 0, 2, 4, 6, 0, 3, 6, 9] else: - matrix_test = [0, 1, 2, 3] + if order == "system" or order == "column": + matrix_test = [0, 2, 1, 3] + else: + matrix_test = [0, 1, 2, 3] elif nqubits == 2: if order == "row": matrix_test = np.arange(dim**2) @@ -157,9 +161,11 @@ def test_vectorization(backend, nqubits, order): matrix_test = backend.cast(matrix_test) dim = 2**nqubits - matrix = np.arange(dim**2).reshape((dim, dim)) - matrix = vectorization(matrix, order, backend=backend) - + if statevector and nqubits == 1: + matrix = np.arange(dim**2) + else: + matrix = np.arange(dim**2).reshape((dim, dim)) + matrix = vectorization(backend.cast(matrix), order, backend=backend) backend.assert_allclose(matrix, matrix_test, atol=PRECISION_TOL) @@ -247,8 +253,10 @@ def chi_superop(pauli_order): @pytest.mark.parametrize("order", ["row", "column"]) def test_to_choi(backend, order): - choi_a0 = to_choi(test_a0, order=order, backend=backend) - choi_a1 = to_choi(test_a1, order=order, backend=backend) + test_a0_ = backend.cast(test_a0) + test_a1_ = backend.cast(test_a1) + choi_a0 = to_choi(test_a0_, order=order, backend=backend) + choi_a1 = to_choi(test_a1_, order=order, backend=backend) choi = choi_a0 + choi_a1 axes = [1, 2] if order == "row" else [0, 3] @@ -261,8 +269,10 @@ def test_to_choi(backend, order): @pytest.mark.parametrize("order", ["row", "column"]) def test_to_liouville(backend, order): - liouville_a0 = to_liouville(test_a0, order=order, backend=backend) - liouville_a1 = to_liouville(test_a1, order=order, backend=backend) + test_a0_ = backend.cast(test_a0) + test_a1_ = backend.cast(test_a1) + liouville_a0 = to_liouville(test_a0_, order=order, backend=backend) + liouville_a1 = to_liouville(test_a1_, order=order, backend=backend) liouville = liouville_a0 + liouville_a1 backend.assert_allclose(liouville, test_superop, atol=PRECISION_TOL) @@ -272,15 +282,17 @@ def test_to_liouville(backend, order): @pytest.mark.parametrize("order", ["row", "column"]) @pytest.mark.parametrize("normalize", [False, True]) def test_to_pauli_liouville(backend, normalize, order, pauli_order): + test_a0_ = backend.cast(test_a0) + test_a1_ = backend.cast(test_a1) pauli_a0 = to_pauli_liouville( - test_a0, + test_a0_, normalize=normalize, order=order, pauli_order=pauli_order, backend=backend, ) pauli_a1 = to_pauli_liouville( - test_a1, + test_a1_, normalize=normalize, order=order, pauli_order=pauli_order, @@ -300,15 +312,17 @@ def test_to_pauli_liouville(backend, normalize, order, pauli_order): @pytest.mark.parametrize("order", ["row", "column"]) @pytest.mark.parametrize("normalize", [False, True]) def test_to_chi(backend, normalize, order, pauli_order): + test_a0_ = backend.cast(test_a0) + test_a1_ = backend.cast(test_a1) chi_a0 = to_chi( - test_a0, + test_a0_, normalize=normalize, order=order, pauli_order=pauli_order, backend=backend, ) chi_a1 = to_chi( - test_a1, + test_a1_, normalize=normalize, order=order, pauli_order=pauli_order, @@ -352,9 +366,9 @@ def test_choi_to_pauli(backend, normalize, order, pauli_order, test_superop): test_superop = backend.cast(test_superop, dtype=test_superop.dtype) axes = [1, 2] if order == "row" else [0, 3] - test_choi = np.swapaxes(np.reshape(test_superop * aux, [2] * 4), *axes).reshape( - [4, 4] - ) + test_choi = backend.np.swapaxes( + backend.np.reshape(test_superop * aux, [2] * 4), *axes + ).reshape([4, 4]) pauli_op = choi_to_pauli( test_choi, normalize, order, pauli_order=pauli_order, backend=backend @@ -401,11 +415,11 @@ def test_choi_to_kraus( state = random_density_matrix(2, backend=backend) - evolution_a0 = a0 @ state @ np.transpose(np.conj(a0)) - evolution_a1 = a1 @ state @ np.transpose(np.conj(a1)) + evolution_a0 = a0 @ state @ backend.np.conj(a0).T + evolution_a1 = a1 @ state @ backend.np.conj(a1).T - test_evolution_a0 = test_a0 @ state @ np.transpose(np.conj(test_a0)) - test_evolution_a1 = test_a1 @ state @ np.transpose(np.conj(test_a1)) + test_evolution_a0 = test_a0 @ state @ backend.np.conj(test_a0).T + test_evolution_a1 = test_a1 @ state @ backend.np.conj(test_a1).T backend.assert_allclose(evolution_a0, test_evolution_a0, atol=2 * PRECISION_TOL) backend.assert_allclose(evolution_a1, test_evolution_a1, atol=2 * PRECISION_TOL) @@ -423,9 +437,9 @@ def test_choi_to_kraus( test_coefficients, ): state = random_density_matrix(2, backend=backend) - evolution = left @ state @ np.transpose(np.conj(right)) + evolution = left @ state @ backend.np.conj(right).T test_evolution = ( - test_coeff**2 * test_left @ state @ np.transpose(np.conj(test_right)) + test_coeff**2 * test_left @ state @ backend.np.conj(test_right).T ) backend.assert_allclose(evolution, test_evolution, atol=2 * PRECISION_TOL) @@ -444,9 +458,9 @@ def test_choi_to_chi(backend, normalize, order, pauli_order, test_superop): test_chi = backend.cast(test_chi, dtype=test_chi.dtype) axes = [1, 2] if order == "row" else [0, 3] - test_choi = np.swapaxes(np.reshape(test_superop * aux, [2] * 4), *axes).reshape( - [4, 4] - ) + test_choi = backend.np.swapaxes( + backend.np.reshape(test_superop * aux, [2] * 4), *axes + ).reshape([4, 4]) chi_matrix = choi_to_chi( test_choi, @@ -506,20 +520,22 @@ def test_choi_to_stinespring( # action of stinespring channel on state + environment stinespring = ( stinespring - @ np.kron(state, np.outer(v_0, v_0)) - @ np.transpose(np.conj(stinespring)) + @ backend.np.kron(state, backend.np.outer(v_0, v_0)) + @ backend.np.conj(stinespring).T ) # partial trace of the environment - stinespring = np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) - stinespring = np.swapaxes(stinespring, 1, 2) + stinespring = backend.np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) + stinespring = backend.np.swapaxes(stinespring, 1, 2) state_final = np.zeros((2**nqubits, 2**nqubits), dtype=complex) state_final = backend.cast(state_final, dtype=state_final.dtype) for alpha in range(2): vector_alpha = np.zeros(2, dtype=complex) vector_alpha[alpha] = 1.0 vector_alpha = backend.cast(vector_alpha, dtype=vector_alpha.dtype) - state_final += np.conj(vector_alpha) @ stinespring @ vector_alpha + state_final = ( + state_final + backend.np.conj(vector_alpha) @ stinespring @ vector_alpha + ) backend.assert_allclose(state_final, state_final_test, atol=PRECISION_TOL) @@ -532,9 +548,9 @@ def test_liouville_to_choi(backend, order, test_superop): choi = liouville_to_choi(test_superop, order=order, backend=backend) axes = [1, 2] if order == "row" else [0, 3] - test_choi = np.reshape(test_superop, [2] * 4) - test_choi = np.swapaxes(test_choi, *axes) - test_choi = np.reshape(test_choi, (4, 4)) + test_choi = backend.np.reshape(test_superop, [2] * 4) + test_choi = backend.np.swapaxes(test_choi, *axes) + test_choi = backend.np.reshape(test_choi, (4, 4)) backend.assert_allclose(choi, test_choi, atol=PRECISION_TOL) @@ -582,11 +598,11 @@ def test_liouville_to_kraus(backend, order, test_a0, test_a1): state = random_density_matrix(2, backend=backend) - evolution_a0 = a0 @ state @ np.transpose(np.conj(a0)) - evolution_a1 = a1 @ state @ np.transpose(np.conj(a1)) + evolution_a0 = a0 @ state @ backend.np.conj(a0).T + evolution_a1 = a1 @ state @ backend.np.conj(a1).T - test_evolution_a0 = test_a0 @ state @ np.transpose(np.conj(test_a0)) - test_evolution_a1 = test_a1 @ state @ np.transpose(np.conj(test_a1)) + test_evolution_a0 = test_a0 @ state @ backend.np.conj(test_a0).T + test_evolution_a1 = test_a1 @ state @ backend.np.conj(test_a1).T backend.assert_allclose(evolution_a0, test_evolution_a0, atol=PRECISION_TOL) backend.assert_allclose(evolution_a1, test_evolution_a1, atol=PRECISION_TOL) @@ -654,20 +670,22 @@ def test_liouville_to_stinespring( # action of stinespring channel on state + environment stinespring = ( stinespring - @ np.kron(state, np.outer(v_0, v_0)) - @ np.transpose(np.conj(stinespring)) + @ backend.np.kron(state, backend.np.outer(v_0, v_0)) + @ backend.np.conj(stinespring).T ) # partial trace of the environment - stinespring = np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) - stinespring = np.swapaxes(stinespring, 1, 2) + stinespring = backend.np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) + stinespring = backend.np.swapaxes(stinespring, 1, 2) state_final = np.zeros((2**nqubits, 2**nqubits), dtype=complex) state_final = backend.cast(state_final, dtype=state_final.dtype) for alpha in range(2): vector_alpha = np.zeros(2, dtype=complex) vector_alpha[alpha] = 1.0 vector_alpha = backend.cast(vector_alpha, dtype=vector_alpha.dtype) - state_final += np.conj(vector_alpha) @ stinespring @ vector_alpha + state_final = ( + state_final + backend.np.conj(vector_alpha) @ stinespring @ vector_alpha + ) backend.assert_allclose(state_final, state_final_test, atol=PRECISION_TOL) @@ -797,8 +815,8 @@ def test_pauli_to_choi(backend, normalize, order, pauli_order, test_superop): ) axes = [1, 2] if order == "row" else [0, 3] - test_choi = np.swapaxes(np.reshape(test_superop, [2] * 4), *axes) - test_choi = np.reshape(test_choi, (4, 4)) + test_choi = backend.np.swapaxes(backend.np.reshape(test_superop, [2] * 4), *axes) + test_choi = backend.np.reshape(test_choi, (4, 4)) backend.assert_allclose(test_choi, choi_super_op, atol=PRECISION_TOL) @@ -833,11 +851,11 @@ def test_pauli_to_kraus(backend, normalize, order, pauli_order, test_a0, test_a1 state = random_density_matrix(2, backend=backend) - evolution_a0 = a0 @ state @ np.transpose(np.conj(a0)) - evolution_a1 = a1 @ state @ np.transpose(np.conj(a1)) + evolution_a0 = a0 @ state @ backend.np.conj(a0).T + evolution_a1 = a1 @ state @ backend.np.conj(a1).T - test_evolution_a0 = test_a0 @ state @ np.transpose(np.conj(test_a0)) - test_evolution_a1 = test_a1 @ state @ np.transpose(np.conj(test_a1)) + test_evolution_a0 = test_a0 @ state @ backend.np.conj(test_a0).T + test_evolution_a1 = test_a1 @ state @ backend.np.conj(test_a1).T backend.assert_allclose(evolution_a0, test_evolution_a0, atol=2 * PRECISION_TOL) backend.assert_allclose(evolution_a1, test_evolution_a1, atol=2 * PRECISION_TOL) @@ -909,20 +927,22 @@ def test_pauli_to_stinespring( # action of stinespring channel on state + environment stinespring = ( stinespring - @ np.kron(state, np.outer(v_0, v_0)) - @ np.transpose(np.conj(stinespring)) + @ backend.np.kron(state, backend.np.outer(v_0, v_0)) + @ backend.np.conj(stinespring).T ) # partial trace of the environment - stinespring = np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) - stinespring = np.swapaxes(stinespring, 1, 2) + stinespring = backend.np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) + stinespring = backend.np.swapaxes(stinespring, 1, 2) state_final = np.zeros((2**nqubits, 2**nqubits), dtype=complex) state_final = backend.cast(state_final, dtype=state_final.dtype) for alpha in range(2): vector_alpha = np.zeros(2, dtype=complex) vector_alpha[alpha] = 1.0 vector_alpha = backend.cast(vector_alpha, dtype=vector_alpha.dtype) - state_final += np.conj(vector_alpha) @ stinespring @ vector_alpha + state_final = ( + state_final + backend.np.conj(vector_alpha) @ stinespring @ vector_alpha + ) backend.assert_allclose(state_final, aux * state_final_test, atol=PRECISION_TOL) @@ -940,7 +960,9 @@ def test_chi_to_choi(backend, normalize, order, pauli_order, test_superop): test_superop = backend.cast(test_superop, dtype=backend.dtype) axes = [1, 2] if order == "row" else [0, 3] - test_choi = np.swapaxes(np.reshape(test_superop, [2] * 4), *axes).reshape([4, 4]) + test_choi = backend.np.swapaxes( + backend.np.reshape(test_superop, [2] * 4), *axes + ).reshape([4, 4]) choi_super_op = chi_to_choi( test_chi / aux, normalize, order, pauli_order, backend=backend @@ -1013,11 +1035,11 @@ def test_chi_to_kraus(backend, normalize, order, pauli_order, test_a0, test_a1): state = random_density_matrix(2, backend=backend) - evolution_a0 = a0 @ state @ np.transpose(np.conj(a0)) - evolution_a1 = a1 @ state @ np.transpose(np.conj(a1)) + evolution_a0 = a0 @ state @ backend.np.conj(a0).T + evolution_a1 = a1 @ state @ backend.np.conj(a1).T - test_evolution_a0 = test_a0 @ state @ np.transpose(np.conj(test_a0)) - test_evolution_a1 = test_a1 @ state @ np.transpose(np.conj(test_a1)) + test_evolution_a0 = test_a0 @ state @ backend.np.conj(test_a0).T + test_evolution_a1 = test_a1 @ state @ backend.np.conj(test_a1).T backend.assert_allclose(evolution_a0, test_evolution_a0, atol=2 * PRECISION_TOL) backend.assert_allclose(evolution_a1, test_evolution_a1, atol=2 * PRECISION_TOL) @@ -1062,20 +1084,22 @@ def test_chi_to_stinespring( # action of stinespring channel on state + environment stinespring = ( stinespring - @ np.kron(state, np.outer(v_0, v_0)) - @ np.transpose(np.conj(stinespring)) + @ backend.np.kron(state, backend.np.outer(v_0, v_0)) + @ backend.np.conj(stinespring).T ) # partial trace of the environment - stinespring = np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) - stinespring = np.swapaxes(stinespring, 1, 2) + stinespring = backend.np.reshape(stinespring, (2**nqubits, 2, 2**nqubits, 2)) + stinespring = backend.np.swapaxes(stinespring, 1, 2) state_final = np.zeros((2**nqubits, 2**nqubits), dtype=complex) state_final = backend.cast(state_final, dtype=state_final.dtype) for alpha in range(2): vector_alpha = np.zeros(2, dtype=complex) vector_alpha[alpha] = 1.0 vector_alpha = backend.cast(vector_alpha, dtype=vector_alpha.dtype) - state_final += np.conj(vector_alpha) @ stinespring @ vector_alpha + state_final = ( + state_final + backend.np.conj(vector_alpha) @ stinespring @ vector_alpha + ) backend.assert_allclose(state_final, aux * state_final_test, atol=PRECISION_TOL) @@ -1272,12 +1296,14 @@ def test_reshuffling(backend, order, test_superop): test_superop = backend.cast(test_superop, dtype=test_superop.dtype) backend.assert_allclose( - np.linalg.norm(reshuffled - backend.cast(test_superop)) < PRECISION_TOL, True + np.linalg.norm(backend.to_numpy(reshuffled) - backend.to_numpy(test_superop)) + < PRECISION_TOL, + True, ) axes = [1, 2] if order == "row" else [0, 3] - test_choi = np.reshape( - np.swapaxes(np.reshape(test_superop, [2] * 4), *axes), [4, 4] + test_choi = backend.np.reshape( + backend.np.swapaxes(backend.np.reshape(test_superop, [2] * 4), *axes), [4, 4] ) reshuffled = _reshuffling(test_choi, order, backend=backend) diff --git a/tests/test_quantum_info_utils.py b/tests/test_quantum_info_utils.py index bc93c230cf..8a375e0f21 100644 --- a/tests/test_quantum_info_utils.py +++ b/tests/test_quantum_info_utils.py @@ -162,9 +162,12 @@ def test_hellinger(backend, validate, kind): prob_q = np.random.rand(10) prob_p /= np.sum(prob_p) prob_q /= np.sum(prob_q) + prob_p = backend.cast(prob_p, dtype=prob_p.dtype) + prob_q = backend.cast(prob_q, dtype=prob_q.dtype) target = float( - backend.calculate_norm(np.sqrt(prob_p) - np.sqrt(prob_q)) / np.sqrt(2) + backend.calculate_norm(backend.np.sqrt(prob_p) - backend.np.sqrt(prob_q)) + / np.sqrt(2) ) prob_p = ( @@ -248,6 +251,6 @@ def test_pqc_integral(backend): pqc_int = pqc_integral(circuit, power_t, samples, backend=backend) - fid = fidelity(pqc_int, pqc_int) + fid = fidelity(pqc_int, pqc_int, backend=backend) backend.assert_allclose(fid, 1.0, atol=PRECISION_TOL) diff --git a/tests/test_transpiler_decompositions.py b/tests/test_transpiler_decompositions.py index 5ce8d97cca..663f73f5ff 100644 --- a/tests/test_transpiler_decompositions.py +++ b/tests/test_transpiler_decompositions.py @@ -12,12 +12,12 @@ def assert_matrices_allclose(gate, natives, backend): target_matrix = gate.matrix(backend) - target_matrix = backend.cast(target_matrix, dtype=target_matrix.dtype) # Remove global phase from target matrix normalisation = np.power( - np.linalg.det(target_matrix), 1 / float(target_matrix.shape[0]), dtype=complex + np.linalg.det(backend.to_numpy(target_matrix)), + 1 / float(target_matrix.shape[0]), + dtype=complex, ) - normalisation = backend.cast(normalisation, dtype=normalisation.dtype) target_unitary = target_matrix / normalisation circuit = Circuit(len(gate.qubits)) @@ -25,15 +25,20 @@ def assert_matrices_allclose(gate, natives, backend): native_matrix = circuit.unitary(backend) # Remove global phase from native matrix normalisation = np.power( - np.linalg.det(native_matrix), 1 / float(native_matrix.shape[0]), dtype=complex + np.linalg.det(backend.to_numpy(native_matrix)), + 1 / float(native_matrix.shape[0]), + dtype=complex, ) - normalisation = backend.cast(normalisation, dtype=normalisation.dtype) native_unitary = native_matrix / normalisation # There can still be phase differences of -1, -1j, 1j c = 0 for phase in [1, -1, 1j, -1j]: - if np.allclose(phase * native_unitary, target_unitary, atol=1e-12): + if np.allclose( + backend.to_numpy(phase * native_unitary), + backend.to_numpy(target_unitary), + atol=1e-12, + ): c = 1 backend.assert_allclose(c, 1) assert_decomposition(circuit, natives) diff --git a/tests/test_transpiler_unitary_decompositions.py b/tests/test_transpiler_unitary_decompositions.py index 6b6c0f4e3e..69e1ffd6d6 100644 --- a/tests/test_transpiler_unitary_decompositions.py +++ b/tests/test_transpiler_unitary_decompositions.py @@ -20,20 +20,21 @@ ) -def bell_unitary(hx, hy, hz): +def bell_unitary(hx, hy, hz, backend): ham = ( - hx * np.kron(matrices.X, matrices.X) - + hy * np.kron(matrices.Y, matrices.Y) - + hz * np.kron(matrices.Z, matrices.Z) + hx * backend.cast(np.kron(matrices.X, matrices.X)) + + hy * backend.cast(np.kron(matrices.Y, matrices.Y)) + + hz * backend.cast(np.kron(matrices.Z, matrices.Z)) ) - return expm(-1j * ham) + return backend.cast(expm(-1j * backend.to_numpy(ham))) def assert_single_qubits(backend, psi, ua, ub): """Assert UA, UB map the maximally entangled basis ``psi`` to the magic basis.""" - uaub = np.kron(ua, ub) + uaub = backend.to_numpy(backend.np.kron(ua, ub)) + psi = backend.to_numpy(psi) for i, j in zip(range(4), [0, 1, 3, 2]): - final_state = np.dot(uaub, psi[:, i]) + final_state = np.matmul(uaub, psi[:, i]) target_state = magic_basis[:, j] fidelity = np.abs(np.dot(np.conj(target_state), final_state)) backend.assert_allclose(fidelity, 1) @@ -65,10 +66,10 @@ def test_eigenbasis_entanglement(backend, seed): """Check that the eigenvectors of UT_U are maximally entangled.""" states, eigvals = calculate_psi(unitary, backend=backend) eigvals = backend.cast(eigvals, dtype=eigvals.dtype) - backend.assert_allclose(np.abs(eigvals), np.ones(4)) - for state in np.transpose(states): + backend.assert_allclose(backend.np.abs(eigvals), np.ones(4)) + for state in states.T: state = backend.partial_trace(state, [1], 2) - backend.assert_allclose(purity(state), 0.5) + backend.assert_allclose(purity(state, backend=backend), 0.5) @pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)]) @@ -93,7 +94,9 @@ def test_u_decomposition(backend, seed): calculate_psi(unitary, backend=backend) else: psi, eigvals = calculate_psi(unitary, backend=backend) - psi_tilde = np.conj(np.sqrt(eigvals)) * np.dot(unitary, psi) + psi_tilde = backend.np.conj(backend.np.sqrt(eigvals)) * backend.np.matmul( + unitary, psi + ) ua_dagger, ub_dagger = calculate_single_qubit_unitaries( psi_tilde, backend=backend ) @@ -133,14 +136,14 @@ def test_calculate_h_vector(backend, seed): _, _, ud, _, _ = magic_decomposition(unitary, backend=backend) ud_diag = to_bell_diagonal(ud, backend=backend) assert ud_diag is not None - hx, hy, hz = calculate_h_vector(ud_diag) - target_matrix = bell_unitary(hx, hy, hz) + hx, hy, hz = calculate_h_vector(ud_diag, backend=backend) + target_matrix = bell_unitary(hx, hy, hz, backend) backend.assert_allclose(ud, target_matrix, atol=PRECISION_TOL) def test_cnot_decomposition(backend): hx, hy, hz = np.random.random(3) - target_matrix = bell_unitary(hx, hy, hz) + target_matrix = bell_unitary(hx, hy, hz, backend) c = Circuit(2) c.add(cnot_decomposition(0, 1, hx, hy, hz, backend)) final_matrix = c.unitary(backend) @@ -149,7 +152,7 @@ def test_cnot_decomposition(backend): def test_cnot_decomposition_light(backend): hx, hy = np.random.random(2) - target_matrix = bell_unitary(hx, hy, 0) + target_matrix = bell_unitary(hx, hy, 0, backend) c = Circuit(2) c.add(cnot_decomposition_light(0, 1, hx, hy, backend)) final_matrix = c.unitary(backend) @@ -195,7 +198,7 @@ def test_two_qubit_decomposition_bell_unitary(backend, hz_zero): hx, hy, hz = (2 * np.random.random(3) - 1) * np.pi if hz_zero: hz = 0 - unitary = backend.cast(bell_unitary(hx, hy, hz)) + unitary = backend.cast(bell_unitary(hx, hy, hz, backend)) c = Circuit(2) c.add(two_qubit_decomposition(0, 1, unitary, backend=backend)) final_matrix = c.unitary(backend)