diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index fec993bdb4..81711c0bbe 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -93,7 +93,7 @@ def identity_density_matrix(self, nqubits, normalize: bool = True): def plus_state(self, nqubits): state = self.np.ones(2**nqubits, dtype=self.dtype) - state /= self.np.sqrt(2**nqubits) + state /= self.np.sqrt(self.cast(2**nqubits)) return state def plus_density_matrix(self, nqubits): @@ -114,7 +114,7 @@ def matrix_parametrized(self, gate): def matrix_fused(self, fgate): rank = len(fgate.target_qubits) - matrix = np.eye(2**rank, dtype=self.dtype) + matrix = np.eye(2**rank, dtype=np.complex128) for gate in fgate.gates: # transfer gate matrix to numpy as it is more efficient for # small tensor calculations @@ -122,7 +122,7 @@ def matrix_fused(self, fgate): gmatrix = self.to_numpy(gate.matrix(self)) # Kronecker product with identity is needed to make the # original matrix have shape (2**rank x 2**rank) - eye = np.eye(2 ** (rank - len(gate.qubits)), dtype=self.dtype) + eye = np.eye(2 ** (rank - len(gate.qubits)), dtype=np.complex128) gmatrix = np.kron(gmatrix, eye) # Transpose the new matrix indices so that it targets the # target qubits of the original gate @@ -137,7 +137,7 @@ def matrix_fused(self, fgate): gmatrix = np.reshape(gmatrix, original_shape) # fuse the individual gate matrix to the total ``FusedGate`` matrix matrix = gmatrix @ matrix - return matrix + return self.cast(matrix) def control_matrix(self, gate): if len(gate.control_qubits) > 1: @@ -530,7 +530,7 @@ def execute_circuit_repeated(self, circuit, nshots, initial_state=None): if circuit.density_matrix: # this implies also it has_collapse assert circuit.has_collapse - final_state = self.np.mean(self.to_numpy(final_states), 0) + final_state = np.mean(self.to_numpy(final_states), 0) if circuit.measurements: qubits = [q for m in circuit.measurements for q in m.target_qubits] final_result = CircuitResult( @@ -741,10 +741,9 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): else: from scipy.linalg import expm return expm(-1j * a * matrix) - else: - expd = self.np.diag(self.np.exp(-1j * a * eigenvalues)) - ud = self.np.transpose(self.np.conj(eigenvectors)) - return self.np.matmul(eigenvectors, self.np.matmul(expd, ud)) + expd = self.np.diag(self.np.exp(-1j * a * eigenvalues)) + ud = self.np.transpose(self.np.conj(eigenvectors)) + return self.np.matmul(eigenvectors, self.np.matmul(expd, ud)) def calculate_expectation_state(self, hamiltonian, state, normalize): statec = self.np.conj(state) diff --git a/src/qibo/backends/pytorch.py b/src/qibo/backends/pytorch.py index 0f3f781046..e34f7d290c 100644 --- a/src/qibo/backends/pytorch.py +++ b/src/qibo/backends/pytorch.py @@ -1,4 +1,3 @@ -import collections from typing import Union import numpy as np @@ -58,6 +57,9 @@ def __init__(self): def set_device(self, device): # pragma: no cover self.device = device + def set_seed(self, seed): + self.np.manual_seed(seed) + def cast( self, x: Union[torch.Tensor, list[torch.Tensor], np.ndarray, list[np.ndarray]], @@ -125,6 +127,8 @@ def issparse(self, x): return super().issparse(x) def to_numpy(self, x): + if isinstance(x, list): + return np.asarray([self.to_numpy(i) for i in x]) if isinstance(x, self.np.Tensor): return x.numpy(force=True) return x @@ -140,10 +144,6 @@ def matrix_parametrized(self, gate): npmatrix = super().matrix_parametrized(gate) return self.np.tensor(npmatrix, dtype=self.dtype) - def matrix_fused(self, gate): - npmatrix = super().matrix_fused(gate) - return self.np.tensor(npmatrix, dtype=self.dtype) - def sample_shots(self, probabilities, nshots): return self.np.multinomial( self.cast(probabilities, dtype="float"), nshots, replacement=True @@ -185,10 +185,10 @@ def calculate_probabilities_density_matrix(self, state, qubits, nqubits): probs = self.np.reshape(probs, len(qubits) * (2,)) return self._order_probabilities(probs, qubits, nqubits).view(-1) - def calculate_frequencies(self, samples): - res, counts = self.np.unique(samples, return_counts=True) - res, counts = res.tolist(), counts.tolist() - return collections.Counter({k: v for k, v in zip(res, counts)}) + # def calculate_frequencies(self, samples): + # res, counts = self.np.unique(samples, return_counts=True) + # res, counts = res.tolist(), counts.tolist() + # return collections.Counter({k: v for k, v in zip(res, counts)}) def update_frequencies(self, frequencies, probabilities, nsamples): frequencies = self.cast(frequencies, dtype="int") @@ -220,7 +220,9 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): return self.np.linalg.matrix_exp( # pylint: disable=not-callable -1j * a * matrix ) - return super().calculate_matrix_exp(a, matrix, eigenvectors, eigenvalues) + expd = self.np.diag(self.np.exp(-1j * a * eigenvalues)) + ud = self.np.transpose(self.np.conj(eigenvectors), 0, 1) + return self.np.matmul(eigenvectors, self.np.matmul(expd, ud)) def calculate_expectation_state(self, hamiltonian, state, normalize): state = self.cast(state) diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py index 71f71f49ed..179895f751 100644 --- a/src/qibo/hamiltonians/hamiltonians.py +++ b/src/qibo/hamiltonians/hamiltonians.py @@ -243,6 +243,7 @@ def __mul__(self, o): ) new_matrix = self.matrix * o r = self.__class__(self.nqubits, new_matrix, backend=self.backend) + o = self.backend.cast(o) if self._eigenvalues is not None: if self.backend.np.real(o) >= 0: # TODO: check for side effects K.qnp r._eigenvalues = o * self._eigenvalues diff --git a/src/qibo/solvers.py b/src/qibo/solvers.py index 8a24ff3c3d..69ab0d3622 100644 --- a/src/qibo/solvers.py +++ b/src/qibo/solvers.py @@ -74,7 +74,7 @@ class Exponential(BaseSolver): def __call__(self, state): propagator = self.current_hamiltonian.exp(self.dt) self.t += self.dt - return (propagator @ state[:, self.backend.np.newaxis])[:, 0] + return (propagator @ state[:, None])[:, 0] class RungeKutta4(BaseSolver): diff --git a/tests/test_backends.py b/tests/test_backends.py index a4714bd378..61ed409f49 100644 --- a/tests/test_backends.py +++ b/tests/test_backends.py @@ -108,7 +108,7 @@ def test_control_matrix_unitary(backend): u = np.random.random((2, 2)) gate = gates.Unitary(u, 0).controlled_by(1) matrix = backend.control_matrix(gate) - target_matrix = np.eye(4, dtype=backend.dtype) + target_matrix = np.eye(4, dtype=np.complex128) target_matrix[2:, 2:] = u backend.assert_allclose(matrix, target_matrix) diff --git a/tests/test_gates_gates.py b/tests/test_gates_gates.py index 09560f4d20..f2d5f1db67 100644 --- a/tests/test_gates_gates.py +++ b/tests/test_gates_gates.py @@ -753,7 +753,7 @@ def test_fswap(backend): [0, 1, 0, 0], [0, 0, 0, -1], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype) target_state = matrix @ initial_state @@ -815,7 +815,7 @@ def test_sycamore(backend): [0, -1j, 0, 0], [0, 0, 0, np.exp(-1j * np.pi / 6)], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype) target_state = matrix @ initial_state @@ -955,7 +955,7 @@ def test_rzx(backend): [0, 0, cos, 1j * sin], [0, 0, 1j * sin, cos], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype) target_state = matrix @ initial_state @@ -996,7 +996,7 @@ def test_rxxyy(backend): [0, -1j * sin, cos, 0], [0, 0, 0, 1], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype) target_state = matrix @ initial_state @@ -1081,7 +1081,7 @@ def test_givens(backend): [0, np.sin(theta), np.cos(theta), 0], [0, 0, 0, 1], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype) @@ -1121,7 +1121,7 @@ def test_rbs(backend): [0, -np.sin(theta), np.cos(theta), 0], [0, 0, 0, 1], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype) @@ -1160,7 +1160,7 @@ def test_ecr(backend): [1, -1j, 0, 0], [-1j, 1, 0, 0], ], - dtype=backend.dtype, + dtype=np.complex128, ) / np.sqrt(2) matrix = backend.cast(matrix, dtype=matrix.dtype) @@ -1224,7 +1224,7 @@ def test_deutsch(backend): [0, 0, 0, 0, 0, 0, 1j * cos, sin], [0, 0, 0, 0, 0, 0, sin, 1j * cos], ], - dtype=backend.dtype, + dtype=np.complex128, ) matrix = backend.cast(matrix, dtype=matrix.dtype)