solved tests and removed duplicate functions

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,15 +114,15 @@ 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
             # explicit to_numpy see https://github.com/qiboteam/qibo/issues/928
             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)

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)

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

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

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)
 

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)