diff --git a/.pylintrc b/.pylintrc
index 7136ff2f74..05d6d81aef 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -10,7 +10,7 @@ fail-under=10
 
 # Add files or directories to the blacklist. They should be base names, not
 # paths.
-ignore=CVS
+ignore=CVS,tests,core,models
 
 # Add files or directories matching the regex patterns to the blacklist. The
 # regex matches against base names, not paths.
diff --git a/src/qibo/backends/__init__.py b/src/qibo/backends/__init__.py
index 8812b8e089..319282fdab 100644
--- a/src/qibo/backends/__init__.py
+++ b/src/qibo/backends/__init__.py
@@ -102,7 +102,9 @@ def set_device(device):
     if device[0] != "/" or len(parts) < 2 or len(parts) > 3:
         raise_error(ValueError, "Device name should follow the pattern: "
                                 "/{device type}:{device number}.")
-    GlobalBackend().set_device(device)
+    backend = GlobalBackend()
+    backend.set_device(device)
+    log.info(f"Using {backend} backend on {backend.device}")
 
 
 def get_threads():
diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py
index fa16f85f35..8a3f54044f 100644
--- a/src/qibo/backends/abstract.py
+++ b/src/qibo/backends/abstract.py
@@ -37,10 +37,6 @@ def apply_gate(self, gate, state, nqubits): # pragma: no cover
     def apply_gate_density_matrix(self, gate, state, nqubits): # pragma: no cover
         raise_error(NotImplementedError)
 
-    @abc.abstractmethod
-    def apply_gate_half_density_matrix(self, gate, state, nqubits): # pragma: no cover
-        raise_error(NotImplementedError)
-
     @abc.abstractmethod
     def execute_circuit(self, circuit, nshots=None): # pragma: no cover
         raise_error(NotImplementedError)
@@ -95,11 +91,15 @@ def cast(self, x, copy=False): # pragma: no cover
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
-    def issparse(self, x):
+    def issparse(self, x): # pragma: no cover
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
-    def to_numpy(self, x):
+    def to_numpy(self, x): # pragma: no cover
+        raise_error(NotImplementedError)
+
+    @abc.abstractmethod
+    def compile(self, func): # pragma: no cover
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
@@ -139,6 +139,10 @@ def apply_gate(self, gate): # pragma: no cover
     def apply_gate_density_matrix(self, gate): # pragma: no cover
         raise_error(NotImplementedError)
 
+    @abc.abstractmethod
+    def apply_gate_half_density_matrix(self, gate, state, nqubits): # pragma: no cover
+        raise_error(NotImplementedError)
+
     @abc.abstractmethod
     def apply_channel(self, gate): # pragma: no cover
         raise_error(NotImplementedError)
@@ -148,14 +152,14 @@ def apply_channel_density_matrix(self, gate): # pragma: no cover
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
-    def collapse_state(self, gate, state, nqubits):
+    def collapse_state(self, state, qubits, shot, nqubits, normalize=True):
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
-    def collapse_density_matrix(self, gate, state, nqubits):
+    def collapse_density_matrix(self, state, qubits, shot, nqubits, normalize=True):
         raise_error(NotImplementedError)
 
-    def execute_circuit(self, circuit, initial_state=None, nshots=None):
+    def execute_circuit(self, circuit, initial_state=None, nshots=None, return_array=False):
         from qibo.gates.special import CallbackGate
 
         if circuit.accelerators and not self.supports_multigpu:
@@ -189,10 +193,12 @@ def execute_circuit(self, circuit, initial_state=None, nshots=None):
                 for gate in circuit.queue:
                     state = gate.apply(self, state, nqubits)
 
-            # TODO: Consider implementing a final state setter in circuits?
-            circuit._final_state = CircuitResult(self, circuit, state, nshots)
-            return circuit._final_state
-        
+            if return_array:
+                return state
+            else:
+                circuit._final_state = CircuitResult(self, circuit, state, nshots)
+                return circuit._final_state
+
         except self.oom_error:
             raise_error(RuntimeError, f"State does not fit in {self.device} memory."
                                        "Please switch the execution device to a "
@@ -290,6 +296,10 @@ def sample_frequencies(self, probabilities, nshots): # pragma: no cover
     def calculate_frequencies(self, samples): # pragma: no cover
         raise_error(NotImplementedError)
 
+    @abc.abstractmethod
+    def update_frequencies(self, frequencies, probabilities, nsamples): # pragma: no cover
+        raise_error(NotImplementedError)
+
     @abc.abstractmethod
     def partial_trace(self, state, qubits, nqubits): # pragma: no cover
         raise_error(NotImplementedError)
diff --git a/src/qibo/backends/matrices.py b/src/qibo/backends/matrices.py
index c3dcefa132..a42342314b 100644
--- a/src/qibo/backends/matrices.py
+++ b/src/qibo/backends/matrices.py
@@ -1,66 +1,78 @@
-import numpy as np
-from functools import cached_property
+import sys
 from qibo.config import raise_error
 
+if sys.version_info.minor >= 8:
+    from functools import cached_property  # pylint: disable=E0611
+else:
+    # Custom ``cached_property`` because it is not available for Python < 3.8
+    from functools import lru_cache
+    def cached_property(func):
+        @property
+        def wrapper(self):
+            return lru_cache()(func)(self)
+        return wrapper
+
 
 class Matrices:
     """Matrix representation of every gate as a numpy array."""
 
     def __init__(self, dtype):
+        import numpy as np
         self.dtype = dtype
+        self.np = np
 
     @cached_property
     def H(self):
-        return np.array([
+        return self.np.array([
             [1, 1], 
             [1, -1]
-        ], dtype=self.dtype) / np.sqrt(2)
+        ], dtype=self.dtype) / self.np.sqrt(2)
 
     @cached_property
     def X(self):
-        return np.array([
+        return self.np.array([
             [0, 1], 
             [1, 0]
         ], dtype=self.dtype)
 
     @cached_property
     def Y(self):
-        return np.array([
+        return self.np.array([
             [0, -1j], 
             [1j, 0]
         ], dtype=self.dtype)
 
     @cached_property
     def Z(self):
-        return np.array([
+        return self.np.array([
             [1, 0], 
             [0, -1]
         ], dtype=self.dtype)
 
     @cached_property
     def S(self):
-        return np.array([
+        return self.np.array([
             [1, 0], 
             [0, 1j]
         ], dtype=self.dtype)
 
     @cached_property
     def SDG(self):
-        return np.conj(self.S)
+        return self.np.conj(self.S)
 
     @cached_property
     def T(self):
-        return np.array([
+        return self.np.array([
             [1, 0],
-            [0, np.exp(1j * np.pi / 4.0)]
+            [0, self.np.exp(1j * self.np.pi / 4.0)]
         ], dtype=self.dtype)
 
     @cached_property
     def TDG(self):
-        return np.conj(self.T)
+        return self.np.conj(self.T)
 
     def I(self, n=2):
-        return np.eye(n, dtype=self.dtype)
+        return self.np.eye(n, dtype=self.dtype)
 
     def Align(self, n=2):
         return self.I(n)
@@ -69,56 +81,56 @@ def M(self):
         raise_error(NotImplementedError)
 
     def RX(self, theta):
-        cos = np.cos(theta / 2.0) + 0j
-        isin = -1j * np.sin(theta / 2.0)
-        return np.array([
+        cos = self.np.cos(theta / 2.0) + 0j
+        isin = -1j * self.np.sin(theta / 2.0)
+        return self.np.array([
             [cos, isin], 
             [isin, cos]
         ], dtype=self.dtype)
 
     def RY(self, theta):
-        cos = np.cos(theta / 2.0) + 0j
-        sin = np.sin(theta / 2.0)
-        return np.array([
+        cos = self.np.cos(theta / 2.0) + 0j
+        sin = self.np.sin(theta / 2.0)
+        return self.np.array([
             [cos, -sin], 
             [sin, cos]
         ], dtype=self.dtype)
 
     def RZ(self, theta):
-        phase = np.exp(0.5j * theta)
-        return np.array([
-            [np.conj(phase), 0], 
+        phase = self.np.exp(0.5j * theta)
+        return self.np.array([
+            [self.np.conj(phase), 0], 
             [0, phase]
         ], dtype=self.dtype)
 
     def U1(self, theta):
-        phase = np.exp(1j * theta)
-        return np.array([
+        phase = self.np.exp(1j * theta)
+        return self.np.array([
             [1, 0], 
             [0, phase]
         ], dtype=self.dtype)
 
     def U2(self, phi, lam):
-        eplus = np.exp(1j * (phi + lam) / 2.0)
-        eminus = np.exp(1j * (phi - lam) / 2.0)
-        return np.array([
-            [np.conj(eplus), - np.conj(eminus)],
+        eplus = self.np.exp(1j * (phi + lam) / 2.0)
+        eminus = self.np.exp(1j * (phi - lam) / 2.0)
+        return self.np.array([
+            [self.np.conj(eplus), - self.np.conj(eminus)],
             [eminus, eplus]
-        ], dtype=self.dtype) / np.sqrt(2)
+        ], dtype=self.dtype) / self.np.sqrt(2)
 
     def U3(self, theta, phi, lam):
-        cost = np.cos(theta / 2)
-        sint = np.sin(theta / 2)
-        eplus = np.exp(1j * (phi + lam) / 2.0)
-        eminus = np.exp(1j * (phi - lam) / 2.0)
-        return np.array([
-            [np.conj(eplus) * cost, - np.conj(eminus) * sint],
+        cost = self.np.cos(theta / 2)
+        sint = self.np.sin(theta / 2)
+        eplus = self.np.exp(1j * (phi + lam) / 2.0)
+        eminus = self.np.exp(1j * (phi - lam) / 2.0)
+        return self.np.array([
+            [self.np.conj(eplus) * cost, - self.np.conj(eminus) * sint],
             [eminus * sint, eplus * cost]
         ], dtype=self.dtype)
 
     @cached_property
     def CNOT(self):
-        return np.array([
+        return self.np.array([
             [1, 0, 0, 0], 
             [0, 1, 0, 0], 
             [0, 0, 0, 1], 
@@ -127,7 +139,7 @@ def CNOT(self):
 
     @cached_property
     def CZ(self):
-        return np.array([
+        return self.np.array([
             [1, 0, 0, 0], 
             [0, 1, 0, 0], 
             [0, 0, 1, 0], 
@@ -135,38 +147,38 @@ def CZ(self):
         ], dtype=self.dtype)
 
     def CRX(self, theta):
-        m = np.eye(4, dtype=self.dtype)
+        m = self.np.eye(4, dtype=self.dtype)
         m[2:, 2:] = self.RX(theta)
         return m
 
     def CRY(self, theta):
-        m = np.eye(4, dtype=self.dtype)
+        m = self.np.eye(4, dtype=self.dtype)
         m[2:, 2:] = self.RY(theta)
         return m
 
     def CRZ(self, theta):
-        m = np.eye(4, dtype=self.dtype)
+        m = self.np.eye(4, dtype=self.dtype)
         m[2:, 2:] = self.RZ(theta)
         return m
 
     def CU1(self, theta):
-        m = np.eye(4, dtype=self.dtype)
+        m = self.np.eye(4, dtype=self.dtype)
         m[2:, 2:] = self.U1(theta)
         return m
 
     def CU2(self, phi, lam):
-        m = np.eye(4, dtype=self.dtype)
+        m = self.np.eye(4, dtype=self.dtype)
         m[2:, 2:] = self.U2(phi, lam)
         return m
 
     def CU3(self, theta, phi, lam):
-        m = np.eye(4, dtype=self.dtype)
+        m = self.np.eye(4, dtype=self.dtype)
         m[2:, 2:] = self.U3(theta, phi, lam)
         return m
 
     @cached_property
     def SWAP(self):
-        return np.array([
+        return self.np.array([
             [1, 0, 0, 0], 
             [0, 0, 1, 0], 
             [0, 1, 0, 0], 
@@ -175,7 +187,7 @@ def SWAP(self):
 
     @cached_property
     def FSWAP(self):
-        return np.array([
+        return self.np.array([
             [1, 0, 0, 0], 
             [0, 0, 1, 0], 
             [0, 1, 0, 0], 
@@ -183,10 +195,10 @@ def FSWAP(self):
         ], dtype=self.dtype)
 
     def fSim(self, theta, phi):
-        cost = np.cos(theta) + 0j
-        isint = -1j * np.sin(theta)
-        phase = np.exp(-1j * phi)
-        return np.array([
+        cost = self.np.cos(theta) + 0j
+        isint = -1j * self.np.sin(theta)
+        phase = self.np.exp(-1j * phi)
+        return self.np.array([
             [1, 0, 0, 0],
             [0, cost, isint, 0],
             [0, isint, cost, 0],
@@ -194,8 +206,8 @@ def fSim(self, theta, phi):
         ], dtype=self.dtype)
 
     def GeneralizedfSim(self, u, phi):
-        phase = np.exp(-1j * phi)
-        return np.array([
+        phase = self.np.exp(-1j * phi)
+        return self.np.array([
             [1, 0, 0, 0],
             [0, u[0, 0], u[0, 1], 0],
             [0, u[1, 0], u[1, 1], 0],
@@ -204,13 +216,19 @@ def GeneralizedfSim(self, u, phi):
 
     @cached_property
     def TOFFOLI(self):
-        m = np.eye(8, dtype=self.dtype)
-        m[-2, -2], m[-2, -1] = 0, 1
-        m[-1, -2], m[-1, -1] = 1, 0
-        return m
+        return self.np.array([
+            [1, 0, 0, 0, 0, 0, 0, 0],
+            [0, 1, 0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0, 0, 0],
+            [0, 0, 0, 1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 1],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+        ])
 
     def Unitary(self, u):
-        return np.array(u, dtype=self.dtype, copy=False)
+        return self.np.array(u, dtype=self.dtype, copy=False)
 
     def VariationalLayer(self, *args):
         raise_error(NotImplementedError)
diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index bac0a4e906..df7f02f718 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -1,5 +1,5 @@
-import collections
 import numpy as np
+import collections
 from qibo.config import raise_error, log
 from qibo.gates import FusedGate
 from qibo.backends import einsum_utils
@@ -11,6 +11,7 @@ class NumpyBackend(Simulator):
 
     def __init__(self):
         super().__init__()
+        self.np = np
         self.name = "numpy"
         self.matrices = Matrices(self.dtype)
         self.tensor_types = np.ndarray
@@ -45,6 +46,9 @@ def to_numpy(self, x):
             return x.toarray()
         return x
 
+    def compile(self, func):
+        return func
+
     def zero_state(self, nqubits):
         state = np.zeros(2 ** nqubits, dtype=self.dtype)
         state[0] = 1
@@ -91,78 +95,78 @@ def control_matrix(self, gate):
         if shape != (2, 2):
             raise_error(ValueError, "Cannot use ``control_unitary`` method on "
                                     "gate matrix of shape {}.".format(shape))
-        zeros = np.zeros((2, 2), dtype=self.dtype)
-        part1 = np.concatenate([np.eye(2, dtype=self.dtype), zeros], axis=0)
-        part2 = np.concatenate([zeros, matrix], axis=0)
-        return np.concatenate([part1, part2], axis=1)
+        zeros = self.np.zeros((2, 2), dtype=self.dtype)
+        part1 = self.np.concatenate([self.np.eye(2, dtype=self.dtype), zeros], axis=0)
+        part2 = self.np.concatenate([zeros, matrix], axis=0)
+        return self.np.concatenate([part1, part2], axis=1)
 
     def apply_gate(self, gate, state, nqubits):
         state = self.cast(state)
-        state = np.reshape(state, nqubits * (2,))
+        state = self.np.reshape(state, nqubits * (2,))
         matrix = gate.asmatrix(self)
         if gate.is_controlled_by:
-            matrix = np.reshape(matrix, 2  * len(gate.target_qubits) * (2,))
+            matrix = self.np.reshape(matrix, 2  * len(gate.target_qubits) * (2,))
             ncontrol = len(gate.control_qubits)
             nactive = nqubits - ncontrol
             order, targets = einsum_utils.control_order(gate, nqubits)
-            state = np.transpose(state, order)
+            state = self.np.transpose(state, order)
             # Apply `einsum` only to the part of the state where all controls
             # are active. This should be `state[-1]`
-            state = np.reshape(state, (2 ** ncontrol,) + nactive * (2,))
+            state = self.np.reshape(state, (2 ** ncontrol,) + nactive * (2,))
             opstring = einsum_utils.apply_gate_string(targets, nactive)
-            updates = np.einsum(opstring, state[-1], matrix)
+            updates = self.np.einsum(opstring, state[-1], matrix)
             # Concatenate the updated part of the state `updates` with the
             # part of of the state that remained unaffected `state[:-1]`.
-            state = np.concatenate([state[:-1], updates[np.newaxis]], axis=0)
-            state = np.reshape(state, nqubits * (2,))
+            state = self.np.concatenate([state[:-1], updates[self.np.newaxis]], axis=0)
+            state = self.np.reshape(state, nqubits * (2,))
             # Put qubit indices back to their proper places
-            state = np.transpose(state, einsum_utils.reverse_order(order))
+            state = self.np.transpose(state, einsum_utils.reverse_order(order))
         else:
-            matrix = np.reshape(matrix, 2  * len(gate.qubits) * (2,))
+            matrix = self.np.reshape(matrix, 2  * len(gate.qubits) * (2,))
             opstring = einsum_utils.apply_gate_string(gate.qubits, nqubits)
-            state = np.einsum(opstring, state, matrix)
-        return np.reshape(state, (2 ** nqubits,))
+            state = self.np.einsum(opstring, state, matrix)
+        return self.np.reshape(state, (2 ** nqubits,))
 
     def apply_gate_density_matrix(self, gate, state, nqubits):
         state = self.cast(state)
-        state = np.reshape(state, 2 * nqubits * (2,))
+        state = self.np.reshape(state, 2 * nqubits * (2,))
         matrix = gate.asmatrix(self)
         if gate.is_controlled_by:
-            matrix = np.reshape(matrix, 2  * len(gate.target_qubits) * (2,))
-            matrixc = np.conj(matrix)
+            matrix = self.np.reshape(matrix, 2  * len(gate.target_qubits) * (2,))
+            matrixc = self.np.conj(matrix)
             ncontrol = len(gate.control_qubits)
             nactive = nqubits - ncontrol
             n = 2 ** ncontrol
 
             order, targets = einsum_utils.control_order_density_matrix(gate, nqubits)
-            state = np.transpose(state, order)
-            state = np.reshape(state, 2 * (n,) + 2 * nactive * (2,))
+            state = self.np.transpose(state, order)
+            state = self.np.reshape(state, 2 * (n,) + 2 * nactive * (2,))
 
             leftc, rightc = einsum_utils.apply_gate_density_matrix_controlled_string(targets, nactive)
             state01 = state[:n - 1, n - 1]
-            state01 = np.einsum(rightc, state01, matrixc)
+            state01 = self.np.einsum(rightc, state01, matrixc)
             state10 = state[n - 1, :n - 1]
-            state10 = np.einsum(leftc, state10, matrix)
+            state10 = self.np.einsum(leftc, state10, matrix)
 
             left, right = einsum_utils.apply_gate_density_matrix_string(targets, nactive)
             state11 = state[n - 1, n - 1]
-            state11 = np.einsum(right, state11, matrixc)
-            state11 = np.einsum(left, state11, matrix)
+            state11 = self.np.einsum(right, state11, matrixc)
+            state11 = self.np.einsum(left, state11, matrix)
 
             state00 = state[range(n - 1)]
             state00 = state00[:, range(n - 1)]
-            state01 = np.concatenate([state00, state01[:, np.newaxis]], axis=1)
-            state10 = np.concatenate([state10, state11[np.newaxis]], axis=0)
-            state = np.concatenate([state01, state10[np.newaxis]], axis=0)
-            state = np.reshape(state, 2 * nqubits * (2,))
-            state = np.transpose(state, einsum_utils.reverse_order(order))
+            state01 = self.np.concatenate([state00, state01[:, self.np.newaxis]], axis=1)
+            state10 = self.np.concatenate([state10, state11[self.np.newaxis]], axis=0)
+            state = self.np.concatenate([state01, state10[self.np.newaxis]], axis=0)
+            state = self.np.reshape(state, 2 * nqubits * (2,))
+            state = self.np.transpose(state, einsum_utils.reverse_order(order))
         else:
-            matrix = np.reshape(matrix, 2 * len(gate.qubits) * (2,))
-            matrixc = np.conj(matrix)
+            matrix = self.np.reshape(matrix, 2 * len(gate.qubits) * (2,))
+            matrixc = self.np.conj(matrix)
             left, right = einsum_utils.apply_gate_density_matrix_string(gate.qubits, nqubits)
-            state = np.einsum(right, state, matrixc)
-            state = np.einsum(left, state, matrix)
-        return np.reshape(state, 2 * (2 ** nqubits,))
+            state = self.np.einsum(right, state, matrixc)
+            state = self.np.einsum(left, state, matrix)
+        return self.np.reshape(state, 2 * (2 ** nqubits,))
 
     def apply_gate_half_density_matrix(self, gate, state, nqubits):
         state = self.cast(state)
@@ -181,7 +185,7 @@ def apply_gate_half_density_matrix(self, gate, state, nqubits):
 
     def apply_channel(self, channel, state, nqubits):
         for coeff, gate in zip(channel.coefficients, channel.gates):
-            if np.random.random() < coeff:
+            if self.np.random.random() < coeff:
                 state = self.apply_gate(gate, state, nqubits)
         return state
 
@@ -195,57 +199,46 @@ def apply_channel_density_matrix(self, channel, state, nqubits):
     def _append_zeros(self, state, qubits, results):
         """Helper method for collapse."""
         for q, r in zip(qubits, results):
-            state = np.expand_dims(state, axis=q)
+            state = self.np.expand_dims(state, axis=q)
             if r:
-                state = np.concatenate([np.zeros_like(state), state], axis=q)
+                state = self.np.concatenate([self.np.zeros_like(state), state], axis=q)
             else:
-                state = np.concatenate([state, np.zeros_like(state)], axis=q)
+                state = self.np.concatenate([state, self.np.zeros_like(state)], axis=q)
         return state
 
-    def collapse_state(self, gate, state, nqubits):
+    def collapse_state(self, state, qubits, shot, nqubits, normalize=True):
         state = self.cast(state)
         shape = state.shape
-        qubits = sorted(gate.target_qubits)
-        # measure and get result
-        probs = self.calculate_probabilities(state, gate.qubits, nqubits)
-        shots = self.sample_shots(probs, 1)
-        shots = self.samples_to_binary(shots, len(qubits))[0]
-        # update the gate's result with the measurement outcome
-        gate.result.backend = self
-        gate.result.append(shots)
-        # collapse state
-        state = np.reshape(state, nqubits * (2,))
+        binshot = self.samples_to_binary(shot, len(qubits))[0]
+        state = self.np.reshape(state, nqubits * (2,))
         order = list(qubits) + [q for q in range(nqubits) if q not in qubits]
-        substate = np.transpose(state, order)[tuple(shots)]
-        norm = np.sqrt(np.sum(np.abs(substate) ** 2))
-        state = substate / norm
-        state = self._append_zeros(state, qubits, shots)
-        return np.reshape(state, shape)
+        state = self.np.transpose(state, order)
+        subshape = (2 ** len(qubits),) + (nqubits - len(qubits)) * (2,)
+        state = self.np.reshape(state, subshape)[int(shot)]
+        if normalize:
+            norm = self.np.sqrt(self.np.sum(self.np.abs(state) ** 2))
+            state = state / norm
+        state = self._append_zeros(state, qubits, binshot)
+        return self.np.reshape(state, shape)
 
-    def collapse_density_matrix(self, gate, state, nqubits):
+    def collapse_density_matrix(self, state, qubits, shot, nqubits, normalize=True):
         state = self.cast(state)
         shape = state.shape
-        qubits = sorted(gate.target_qubits)
-        # measure and get result
-        probs = self.calculate_probabilities_density_matrix(state, gate.qubits, nqubits)
-        shots = self.sample_shots(probs, 1)
-        shots = list(self.samples_to_binary(shots, len(qubits))[0])
-        # update the gate's result with the measurement outcome
-        gate.result.backend = self
-        gate.result.append(shots)
-        # collapse state
+        binshot = list(self.samples_to_binary(shot, len(qubits))[0])
         order = list(qubits) + [q + nqubits for q in qubits]
         order.extend(q for q in range(nqubits) if q not in qubits)
         order.extend(q + nqubits for q in range(nqubits) if q not in qubits)
-        shots = 2 * shots
-        state = np.reshape(state, 2 * nqubits * (2,))
-        substate = np.transpose(state, order)[tuple(shots)]
-        n = 2 ** (len(substate.shape) // 2)
-        norm = np.trace(np.reshape(substate, (n, n)))
-        state = substate / norm
+        state = self.np.reshape(state, 2 * nqubits * (2,))
+        state = self.np.transpose(state, order)
+        subshape = 2 * (2 ** len(qubits),) + 2 * (nqubits - len(qubits)) * (2,)
+        state = self.np.reshape(state, subshape)[int(shot), int(shot)]
+        n = 2 ** (len(state.shape) // 2)
+        if normalize:
+            norm = self.np.trace(self.np.reshape(state, (n, n)))
+            state = state / norm
         qubits = qubits + [q + nqubits for q in qubits]
-        state = self._append_zeros(state, qubits, shots)
-        return np.reshape(state, shape)
+        state = self._append_zeros(state, qubits, 2 * binshot)
+        return self.np.reshape(state, shape)
 
     def reset_error_density_matrix(self, gate, state, nqubits):
         from qibo.gates import X
@@ -254,13 +247,13 @@ def reset_error_density_matrix(self, gate, state, nqubits):
         q = gate.target_qubits[0]
         p0, p1 = gate.coefficients[:2]
         trace = self.partial_trace_density_matrix(state, (q,), nqubits)
-        trace = np.reshape(trace, 2 * (nqubits - 1) * (2,))
+        trace = self.np.reshape(trace, 2 * (nqubits - 1) * (2,))
         zero = self.zero_density_matrix(1)
-        zero = np.tensordot(trace, zero, axes=0)
+        zero = self.np.tensordot(trace, zero, axes=0)
         order = list(range(2 * nqubits - 2))
         order.insert(q, 2 * nqubits - 2)
         order.insert(q + nqubits, 2 * nqubits - 1)
-        zero = np.reshape(np.transpose(zero, order), shape)
+        zero = self.np.reshape(self.np.transpose(zero, order), shape)
         state = (1 - p0 - p1) * state + p0 * zero
         return state + p1 * self.apply_gate_density_matrix(X(q), zero, nqubits)
 
@@ -268,7 +261,7 @@ def thermal_error_density_matrix(self, gate, state, nqubits):
         state = self.cast(state)
         shape = state.shape
         state = self.apply_gate(gate, state.ravel(), 2 * nqubits)
-        return np.reshape(state, shape)
+        return self.np.reshape(state, shape)
 
     def calculate_symbolic(self, state, nqubits, decimals=5, cutoff=1e-10, max_terms=20):
         state = self.to_numpy(state)
@@ -298,8 +291,7 @@ def calculate_symbolic_density_matrix(self, state, nqubits, decimals=5, cutoff=1
                 return terms
         return terms
 
-    @staticmethod
-    def _order_probabilities(probs, qubits, nqubits):
+    def _order_probabilities(self, probs, qubits, nqubits):
         """Arrange probabilities according to the given ``qubits`` ordering."""
         unmeasured, reduced = [], {}
         for i in range(nqubits):
@@ -307,117 +299,117 @@ def _order_probabilities(probs, qubits, nqubits):
                 reduced[i] = i - len(unmeasured)
             else:
                 unmeasured.append(i)
-        return np.transpose(probs, [reduced.get(i) for i in qubits])
+        return self.np.transpose(probs, [reduced.get(i) for i in qubits])
 
     def calculate_probabilities(self, state, qubits, nqubits):
-        rtype = state.real.dtype
+        rtype = self.np.real(state).dtype
         unmeasured_qubits = tuple(i for i in range(nqubits) if i not in qubits)
-        state = np.reshape(np.abs(state) ** 2, nqubits * (2,))
-        probs = np.sum(state.astype(rtype), axis=unmeasured_qubits)
+        state = self.np.reshape(self.np.abs(state) ** 2, nqubits * (2,))
+        probs = self.np.sum(state.astype(rtype), axis=unmeasured_qubits)
         return self._order_probabilities(probs, qubits, nqubits).ravel()
 
     def calculate_probabilities_density_matrix(self, state, qubits, nqubits):
-        rtype = state.real.dtype
+        rtype = self.np.real(state).dtype
         order = tuple(sorted(qubits))
         order += tuple(i for i in range(nqubits) if i not in qubits)
         order = order + tuple(i + nqubits for i in order)
         shape = 2 * (2 ** len(qubits), 2 ** (nqubits - len(qubits)))
-        state = np.reshape(state, 2 * nqubits * (2,))
-        state = np.reshape(np.transpose(state, order), shape)
-        probs = np.einsum("abab->a", state).astype(rtype)
-        probs = np.reshape(probs, len(qubits) * (2,))
+        state = self.np.reshape(state, 2 * nqubits * (2,))
+        state = self.np.reshape(self.np.transpose(state, order), shape)
+        probs = self.np.einsum("abab->a", state).astype(rtype)
+        probs = self.np.reshape(probs, len(qubits) * (2,))
         return self._order_probabilities(probs, qubits, nqubits).ravel()
 
     def set_seed(self, seed):
-        np.random.seed(seed)
+        self.np.random.seed(seed)
 
     def sample_shots(self, probabilities, nshots):
-        return np.random.choice(range(len(probabilities)), size=nshots, p=probabilities)
+        return self.np.random.choice(range(len(probabilities)), size=nshots, p=probabilities)
 
     def aggregate_shots(self, shots):
-        return np.array(shots, dtype=shots[0].dtype)
+        return self.np.array(shots, dtype=shots[0].dtype)
 
     def samples_to_binary(self, samples, nqubits):
-        qrange = np.arange(nqubits - 1, -1, -1, dtype="int32")
-        return np.mod(np.right_shift(samples[:, np.newaxis], qrange), 2)
+        qrange = self.np.arange(nqubits - 1, -1, -1, dtype="int32")
+        return self.np.mod(self.np.right_shift(samples[:, self.np.newaxis], qrange), 2)
 
     def samples_to_decimal(self, samples, nqubits):
-        qrange = np.arange(nqubits - 1, -1, -1, dtype="int32")
-        qrange = (2 ** qrange)[:, np.newaxis]
-        return np.matmul(samples, qrange)[:, 0]
+        qrange = self.np.arange(nqubits - 1, -1, -1, dtype="int32")
+        qrange = (2 ** qrange)[:, self.np.newaxis]
+        return self.np.matmul(samples, qrange)[:, 0]
+
+    def calculate_frequencies(self, samples):
+        res, counts = self.np.unique(samples, return_counts=True)
+        res, counts = self.np.array(res), self.np.array(counts)
+        return collections.Counter({k: v for k, v in zip(res, counts)})
+
+    def update_frequencies(self, frequencies, probabilities, nsamples):
+        samples = self.sample_shots(probabilities, nsamples)
+        res, counts = self.np.unique(samples, return_counts=True)
+        frequencies[res] += counts
+        return frequencies
 
     def sample_frequencies(self, probabilities, nshots):
         from qibo.config import SHOT_BATCH_SIZE
-        nprobs = probabilities / np.sum(probabilities)
-        def update_frequencies(nsamples, frequencies):
-            samples = np.random.choice(range(len(nprobs)), size=nsamples, p=nprobs)
-            res, counts = np.unique(samples, return_counts=True)
-            frequencies[res] += counts
-            return frequencies
-
-        frequencies = np.zeros(len(nprobs), dtype="int64")
+        nprobs = probabilities / self.np.sum(probabilities)
+        frequencies = self.np.zeros(len(nprobs), dtype="int64")
         for _ in range(nshots // SHOT_BATCH_SIZE):
-            frequencies = update_frequencies(SHOT_BATCH_SIZE, frequencies)
-        frequencies = update_frequencies(nshots % SHOT_BATCH_SIZE, frequencies)
+            frequencies = self.update_frequencies(frequencies, nprobs, SHOT_BATCH_SIZE)
+        frequencies = self.update_frequencies(frequencies, nprobs, nshots % SHOT_BATCH_SIZE)
         return collections.Counter({i: f for i, f in enumerate(frequencies) if f > 0})
 
-    def calculate_frequencies(self, samples):
-        res, counts = np.unique(samples, return_counts=True)
-        res, counts = np.array(res), np.array(counts)
-        return collections.Counter({k: v for k, v in zip(res, counts)})
-
     def apply_bitflips(self, noiseless_samples, bitflip_probabilities):
-        fprobs = np.array(bitflip_probabilities, dtype="float64")
-        sprobs = np.random.random(noiseless_samples.shape)
-        flip0 = np.array(sprobs < fprobs[0], dtype=noiseless_samples.dtype)
-        flip1 = np.array(sprobs < fprobs[1], dtype=noiseless_samples.dtype)
+        fprobs = self.np.array(bitflip_probabilities, dtype="float64")
+        sprobs = self.np.random.random(noiseless_samples.shape)
+        flip0 = self.np.array(sprobs < fprobs[0], dtype=noiseless_samples.dtype)
+        flip1 = self.np.array(sprobs < fprobs[1], dtype=noiseless_samples.dtype)
         noisy_samples = noiseless_samples + (1 - noiseless_samples) * flip0
         noisy_samples = noisy_samples - noiseless_samples * flip1
         return noisy_samples
 
     def partial_trace(self, state, qubits, nqubits):
         state = self.cast(state)
-        state = np.reshape(state, nqubits * (2,))
+        state = self.np.reshape(state, nqubits * (2,))
         axes = 2 * [list(qubits)]
-        rho = np.tensordot(state, np.conj(state), axes=axes)
+        rho = self.np.tensordot(state, self.np.conj(state), axes=axes)
         shape = 2 * (2 ** (nqubits - len(qubits)),)
-        return np.reshape(rho, shape)
+        return self.np.reshape(rho, shape)
 
     def partial_trace_density_matrix(self, state, qubits, nqubits):
         state = self.cast(state)
-        state = np.reshape(state, 2 * nqubits * (2,))
+        state = self.np.reshape(state, 2 * nqubits * (2,))
 
         order = tuple(sorted(qubits))
         order += tuple(i for i in range(nqubits) if i not in qubits)
         order += tuple(i + nqubits for i in order)
         shape = 2 * (2 ** len(qubits), 2 ** (nqubits - len(qubits)))
         
-        state = np.transpose(state, order)
-        state = np.reshape(state, shape)
-        return np.einsum("abac->bc", state)
+        state = self.np.transpose(state, order)
+        state = self.np.reshape(state, shape)
+        return self.np.einsum("abac->bc", state)
 
     def entanglement_entropy(self, rho):
         from qibo.config import EIGVAL_CUTOFF
         # Diagonalize
-        eigvals = np.linalg.eigvalsh(rho).real
+        eigvals = self.np.linalg.eigvalsh(rho).real
         # Treating zero and negative eigenvalues
         masked_eigvals = eigvals[eigvals > EIGVAL_CUTOFF]
-        spectrum = -1 * np.log(masked_eigvals)
-        entropy = np.sum(masked_eigvals * spectrum) / np.log(2.0)
+        spectrum = -1 * self.np.log(masked_eigvals)
+        entropy = self.np.sum(masked_eigvals * spectrum) / self.np.log(2.0)
         return entropy, spectrum
 
     def calculate_norm(self, state):
         state = self.cast(state)
-        return np.sqrt(np.sum(np.abs(state) ** 2))
+        return self.np.sqrt(self.np.sum(self.np.abs(state) ** 2))
 
     def calculate_norm_density_matrix(self, state):
         state = self.cast(state)
-        return np.trace(state)
+        return self.np.trace(state)
 
     def calculate_overlap(self, state1, state2):
         state1 = self.cast(state1)
         state2 = self.cast(state2)
-        return np.abs(np.sum(np.conj(state1) * state2))
+        return self.np.abs(self.np.sum(self.np.conj(state1) * state2))
 
     def calculate_overlap_density_matrix(self, state1, state2):
         raise_error(NotImplementedError)
@@ -445,23 +437,23 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
                 from scipy.linalg import expm
             return expm(-1j * a * matrix)
         else:
-            expd = np.diag(np.exp(-1j * a * eigenvalues))
-            ud = np.transpose(np.conj(eigenvectors))
-            return np.matmul(eigenvectors, 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, matrix, state, normalize):
-        statec = np.conj(state)
+        statec = self.np.conj(state)
         hstate = matrix @ state
-        ev = np.real(np.sum(statec * hstate))
+        ev = self.np.real(self.np.sum(statec * hstate))
         if normalize:
-            norm = np.sum(np.square(np.abs(state)))
+            norm = self.np.sum(self.np.square(self.np.abs(state)))
             ev = ev / norm
         return ev
 
     def calculate_expectation_density_matrix(self, matrix, state, normalize):
-        ev = np.real(np.trace(matrix @ state))
+        ev = self.np.real(self.np.trace(matrix @ state))
         if normalize:
-            norm = np.real(np.trace(state))
+            norm = self.np.real(self.np.trace(state))
             ev = ev / norm
         return ev
 
diff --git a/src/qibo/backends/tensorflow.py b/src/qibo/backends/tensorflow.py
index 6e9d4e4453..f43860e713 100644
--- a/src/qibo/backends/tensorflow.py
+++ b/src/qibo/backends/tensorflow.py
@@ -1,8 +1,145 @@
 import os
+import collections
 import numpy as np
-from qibo.backends import einsum_utils
 from qibo.backends.numpy import NumpyBackend
-from qibo.config import raise_error, TF_LOG_LEVEL
+from qibo.backends.matrices import Matrices
+from qibo.config import log, raise_error, TF_LOG_LEVEL
+
+
+class TensorflowMatrices(Matrices):
+    # Redefine parametrized gate matrices for backpropagation to work
+
+    def __init__(self, dtype):
+        super().__init__(dtype)
+        import tensorflow as tf
+        import tensorflow.experimental.numpy as tnp  # pylint: disable=E0401
+        self.tf = tf
+        self.np = tnp
+
+    def RX(self, theta):
+        cos = self.np.cos(theta / 2.0) + 0j
+        isin = -1j * self.np.sin(theta / 2.0)
+        return self.tf.cast([
+            [cos, isin], 
+            [isin, cos]
+        ], dtype=self.dtype)
+
+    def RY(self, theta):
+        cos = self.np.cos(theta / 2.0) + 0j
+        sin = self.np.sin(theta / 2.0) + 0j
+        return self.tf.cast([
+            [cos, -sin], 
+            [sin, cos]
+        ], dtype=self.dtype)
+
+    def RZ(self, theta):
+        phase = self.np.exp(0.5j * theta)
+        return self.tf.cast([
+            [self.np.conj(phase), 0], 
+            [0, phase]
+        ], dtype=self.dtype)
+
+    def U1(self, theta):
+        phase = self.np.exp(1j * theta)
+        return self.tf.cast([
+            [1, 0], 
+            [0, phase]
+        ], dtype=self.dtype)
+
+    def U2(self, phi, lam):
+        eplus = self.np.exp(1j * (phi + lam) / 2.0)
+        eminus = self.np.exp(1j * (phi - lam) / 2.0)
+        return self.tf.cast([
+            [self.np.conj(eplus), - self.np.conj(eminus)],
+            [eminus, eplus]
+        ], dtype=self.dtype) / self.np.sqrt(2)
+
+    def U3(self, theta, phi, lam):
+        cost = self.np.cos(theta / 2)
+        sint = self.np.sin(theta / 2)
+        eplus = self.np.exp(1j * (phi + lam) / 2.0)
+        eminus = self.np.exp(1j * (phi - lam) / 2.0)
+        return self.tf.cast([
+            [self.np.conj(eplus) * cost, - self.np.conj(eminus) * sint],
+            [eminus * sint, eplus * cost]
+        ], dtype=self.dtype)
+
+    def CRX(self, theta):
+        r = self.RX(theta)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, r[0, 0], r[0, 1]],
+            [0, 0, r[1, 0], r[1, 1]],
+        ], dtype=self.dtype)
+
+    def CRY(self, theta):
+        r = self.RY(theta)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, r[0, 0], r[0, 1]],
+            [0, 0, r[1, 0], r[1, 1]],
+        ], dtype=self.dtype)
+
+    def CRZ(self, theta):
+        r = self.RZ(theta)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, r[0, 0], r[0, 1]],
+            [0, 0, r[1, 0], r[1, 1]],
+        ], dtype=self.dtype)
+
+    def CU1(self, theta):
+        r = self.U1(theta)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, r[0, 0], r[0, 1]],
+            [0, 0, r[1, 0], r[1, 1]],
+        ], dtype=self.dtype)
+
+    def CU2(self, phi, lam):
+        r = self.U2(phi, lam)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, r[0, 0], r[0, 1]],
+            [0, 0, r[1, 0], r[1, 1]],
+        ], dtype=self.dtype)
+
+    def CU3(self, theta, phi, lam):
+        r = self.U3(theta, phi, lam)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, r[0, 0], r[0, 1]],
+            [0, 0, r[1, 0], r[1, 1]],
+        ], dtype=self.dtype)
+
+    def fSim(self, theta, phi):
+        cost = self.np.cos(theta) + 0j
+        isint = -1j * self.np.sin(theta)
+        phase = self.np.exp(-1j * phi)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, cost, isint, 0],
+            [0, isint, cost, 0],
+            [0, 0, 0, phase],
+        ], dtype=self.dtype)
+
+    def GeneralizedfSim(self, u, phi):
+        phase = self.np.exp(-1j * phi)
+        return self.tf.cast([
+            [1, 0, 0, 0],
+            [0, u[0, 0], u[0, 1], 0],
+            [0, u[1, 0], u[1, 1], 0],
+            [0, 0, 0, phase],
+        ], dtype=self.dtype)
+
+    def Unitary(self, u):
+        return self.tf.cast(u, dtype=self.dtype)
 
 
 class TensorflowBackend(NumpyBackend):
@@ -12,51 +149,171 @@ def __init__(self):
         self.name = "tensorflow"
         os.environ["TF_CPP_MIN_LOG_LEVEL"] = str(TF_LOG_LEVEL)
         import tensorflow as tf
+        import tensorflow.experimental.numpy as tnp  # pylint: disable=E0401
+        tnp.experimental_enable_numpy_behavior()
         self.tf = tf
-    
+        self.np = tnp
+        self.matrices = TensorflowMatrices(self.dtype)
+        
+        from tensorflow.python.framework import errors_impl  # pylint: disable=E0611
+        self.oom_error = errors_impl.ResourceExhaustedError
+
+        cpu_devices = tf.config.list_logical_devices("CPU")
+        gpu_devices = tf.config.list_logical_devices("GPU")
+        if gpu_devices: # pragma: no cover
+            # CI does not use GPUs
+            self.device = gpu_devices[0].name
+        elif cpu_devices:
+            self.device = cpu_devices[0].name
+
+        import psutil
+        self.nthreads = psutil.cpu_count(logical=True)
+
+        self.tensor_types = (np.ndarray, tf.Tensor, tf.Variable)
+
     def set_device(self, device):
-        # TODO: Implement this
-        raise_error(NotImplementedError)
+        self.device = device
 
     def set_threads(self, nthreads):
-        # TODO: Implement this
-        raise_error(NotImplementedError)
+        log.warning("`set_threads` is not supported by the tensorflow "
+                    "backend. Please use tensorflow's thread setters: "
+                    "`tf.config.threading.set_inter_op_parallelism_threads` "
+                    "or `tf.config.threading.set_intra_op_parallelism_threads` "
+                    "to switch the number of threads.")
 
-    def asmatrix(self, gate):
-        npmatrix = super().asmatrix(gate)
-        return self.tf.cast(npmatrix, dtype=self.dtype)
+    def cast(self, x, dtype=None, copy=False):
+        if dtype is None:
+            dtype = self.dtype
+        x = self.tf.cast(x, dtype=dtype)
+        if copy:
+            return self.tf.identity(x)
+        return x
 
-    def apply_gate(self, gate, state, nqubits):
-        # TODO: Implement density matrices (most likely in another method)
-        state = self.tf.reshape(state, nqubits * (2,))
-        matrix = self.tf.reshape(self.asmatrix(gate), 2  * len(gate.qubits) * (2,))
-        if gate.is_controlled_by:
-            ncontrol = len(gate.control_qubits)
-            nactive = nqubits - ncontrol
-            order, targets = einsum_utils.control_order(gate, nqubits)
-            state = self.tf.transpose(state, order)
-            # Apply `einsum` only to the part of the state where all controls
-            # are active. This should be `state[-1]`
-            state = self.tf.reshape(state, (2 ** ncontrol,) + nactive * (2,))
-            opstring = einsum_utils.apply_gate_string(targets, nactive)
-            updates = self.tf.einsum(opstring, state[-1], matrix)
-            # Concatenate the updated part of the state `updates` with the
-            # part of of the state that remained unaffected `state[:-1]`.
-            state = self.tf.concatenate([state[:-1], updates[self.tf.newaxis]], axis=0)
-            state = self.tf.reshape(state, nqubits * (2,))
-            # Put qubit indices back to their proper places
-            reverse_order = len(order) * [0]
-            for i, r in enumerate(order):
-                reverse_order[r] = i
-            state = self.tf.transpose(state, reverse_order)
-        else:
-            state = self.tf.einsum(opstring, state, matrix)
-        return self.tf.reshape(state, (2 ** nqubits,))
+    def issparse(self, x):
+        return isinstance(x, self.tf.sparse.SparseTensor)
+
+    def to_numpy(self, x):
+        return np.array(x)
+
+    def compile(self, func):
+        return self.tf.function(func)
 
     def zero_state(self, nqubits):
-        """Generate |000...0> state as an array."""
         idx = self.tf.constant([[0]], dtype="int32")
         state = self.tf.zeros((2 ** nqubits,), dtype=self.dtype)
         update = self.tf.constant([1], dtype=self.dtype)
         state = self.tf.tensor_scatter_nd_update(state, idx, update)
         return state
+
+    def zero_density_matrix(self, nqubits):
+        idx = self.tf.constant([[0, 0]], dtype="int32")
+        state = self.tf.zeros(2 * (2 ** nqubits,), dtype=self.dtype)
+        update = self.tf.constant([1], dtype=self.dtype)
+        state = self.tf.tensor_scatter_nd_update(state, idx, update)
+        return state
+
+    def asmatrix(self, gate):
+        npmatrix = super().asmatrix(gate)
+        return self.tf.cast(npmatrix, dtype=self.dtype)
+
+    def asmatrix_parametrized(self, gate):
+        npmatrix = super().asmatrix_parametrized(gate)
+        return self.tf.cast(npmatrix, dtype=self.dtype)
+
+    def asmatrix_fused(self, gate):
+        npmatrix = super().asmatrix_fused(gate)
+        return self.tf.cast(npmatrix, dtype=self.dtype)
+
+    def execute_circuit(self, circuit, initial_state=None, nshots=None, return_array=False):
+        with self.tf.device(self.device):
+            return super().execute_circuit(circuit, initial_state, nshots, return_array)
+
+    def execute_circuit_repeated(self, circuit, initial_state=None, nshots=None):
+        with self.tf.device(self.device):
+            return super().execute_circuit_repeated(circuit, initial_state, nshots)
+
+    def sample_shots(self, probabilities, nshots):
+        # redefining this because ``tnp.random.choice`` is not available
+        logits = self.tf.math.log(probabilities)[self.tf.newaxis]
+        samples = self.tf.random.categorical(logits, nshots)[0]
+        return samples
+
+    def samples_to_binary(self, samples, nqubits):
+        # redefining this because ``tnp.right_shift`` is not available
+        qrange = self.np.arange(nqubits - 1, -1, -1, dtype="int64")
+        samples = self.tf.bitwise.right_shift(samples[:, self.np.newaxis], qrange)
+        return self.tf.math.mod(samples, 2)
+
+    def calculate_frequencies(self, samples):
+        # redefining this because ``tnp.unique`` is not available
+        res, _, counts = self.tf.unique_with_counts(samples, out_idx="int64")
+        res, counts = self.np.array(res), self.np.array(counts)
+        return collections.Counter({int(k): int(v) for k, v in zip(res, counts)})
+
+    def update_frequencies(self, frequencies, probabilities, nsamples):
+        # redefining this because ``tnp.unique`` and tensor update is not available
+        samples = self.sample_shots(probabilities, nsamples)
+        res, _, counts = self.tf.unique_with_counts(samples, out_idx="int64")
+        frequencies = self.tf.tensor_scatter_nd_add(frequencies, res[:, self.tf.newaxis], counts)
+        return frequencies
+
+    def entanglement_entropy(self, rho):
+        # redefining this because ``tnp.linalg`` is not available
+        from qibo.config import EIGVAL_CUTOFF
+        # Diagonalize
+        eigvals = self.np.real(self.tf.linalg.eigvalsh(rho))
+        # Treating zero and negative eigenvalues
+        masked_eigvals = eigvals[eigvals > EIGVAL_CUTOFF]
+        spectrum = -1 * self.np.log(masked_eigvals)
+        entropy = self.np.sum(masked_eigvals * spectrum) / self.np.log(2.0)
+        return entropy, spectrum
+
+    def calculate_eigenvalues(self, matrix, k=6):
+        return self.tf.linalg.eigvalsh(matrix)
+
+    def calculate_eigenvectors(self, matrix, k=6):
+        return self.tf.linalg.eigh(matrix)
+
+    def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
+        if eigenvectors is None or self.issparse(matrix):
+            return self.tf.linalg.expm(-1j * a * matrix)
+        else:
+            return super().calculate_matrix_exp(a, matrix, eigenvectors, eigenvalues)
+
+    def calculate_matrix_product(self, hamiltonian, o):
+        if isinstance(o, hamiltonian.__class__):
+            new_matrix = self.np.dot(hamiltonian.matrix, o.matrix)
+            return hamiltonian.__class__(hamiltonian.nqubits, new_matrix, backend=self)
+
+        if isinstance(o, self.tensor_types):
+            rank = len(tuple(o.shape))
+            if rank == 1: # vector
+                return self.np.matmul(hamiltonian.matrix, o[:, self.np.newaxis])[:, 0]
+            elif rank == 2: # matrix
+                return self.np.matmul(hamiltonian.matrix, o)
+            else:
+                raise_error(ValueError, "Cannot multiply Hamiltonian with "
+                                        "rank-{} tensor.".format(rank))
+
+        raise_error(NotImplementedError, "Hamiltonian matmul to {} not "
+                                         "implemented.".format(type(o)))
+
+    def test_regressions(self, name):
+        if name == "test_measurementresult_apply_bitflips":
+            return [
+                [4, 0, 0, 1, 0, 2, 2, 4, 4, 0],
+                [4, 0, 0, 1, 0, 2, 2, 4, 4, 0],
+                [4, 0, 0, 1, 0, 0, 0, 4, 4, 0],
+                [4, 0, 0, 0, 0, 0, 0, 4, 4, 0]
+            ]
+        elif name == "test_probabilistic_measurement": 
+            return {0: 271, 1: 239, 2: 242, 3: 248}
+        elif name == "test_unbalanced_probabilistic_measurement": 
+            return {0: 168, 1: 188, 2: 154, 3: 490}
+        elif name == "test_post_measurement_bitflips_on_circuit": 
+            return [
+                {5: 30}, {5: 16, 7: 10, 6: 2, 3: 1, 4: 1},
+                {3: 6, 5: 6, 7: 5, 2: 4, 4: 3, 0: 2, 1: 2, 6: 2}
+            ]
+        else:
+            return None
diff --git a/src/qibo/gates/measurements.py b/src/qibo/gates/measurements.py
index ea47da365c..b310e2e2ed 100644
--- a/src/qibo/gates/measurements.py
+++ b/src/qibo/gates/measurements.py
@@ -214,7 +214,25 @@ def matrix(self):
         raise_error(NotImplementedError, "Measurement gates do not have matrix representation.")
 
     def apply(self, backend, state, nqubits):
-        return backend.collapse_state(self, state, nqubits)
+        qubits = sorted(self.target_qubits)
+        # measure and get result
+        probs = backend.calculate_probabilities(state, qubits, nqubits)
+        shot = backend.sample_shots(probs, 1)
+        # update the gate's result with the measurement outcome
+        binshot = backend.samples_to_binary(shot, len(qubits))[0]
+        self.result.backend = backend
+        self.result.append(binshot)
+        # collapse state
+        return backend.collapse_state(state, qubits, shot, nqubits)
 
     def apply_density_matrix(self, backend, state, nqubits):
-        return backend.collapse_density_matrix(self, state, nqubits)
\ No newline at end of file
+        qubits = sorted(self.target_qubits)
+        # measure and get result
+        probs = backend.calculate_probabilities_density_matrix(state, qubits, nqubits)
+        shot = backend.sample_shots(probs, 1)
+        binshot = backend.samples_to_binary(shot, len(qubits))[0]
+        # update the gate's result with the measurement outcome
+        self.result.backend = backend
+        self.result.append(binshot)
+        # collapse state
+        return backend.collapse_density_matrix(state, qubits, shot, nqubits)
diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py
index f314120a2e..6ba82577a1 100644
--- a/src/qibo/hamiltonians/hamiltonians.py
+++ b/src/qibo/hamiltonians/hamiltonians.py
@@ -28,7 +28,6 @@ def __init__(self, nqubits, matrix=None, backend=None):
         matrix = self.backend.cast(matrix)
 
         super().__init__()
-
         self.nqubits = nqubits
         self.matrix = matrix
         self._eigenvalues = None
@@ -163,12 +162,12 @@ def __mul__(self, o):
         new_matrix = self.matrix * o
         r = self.__class__(self.nqubits, new_matrix, backend=self.backend)
         if self._eigenvalues is not None:
-            if self.backend.cast(o).real >= 0: # TODO: check for side effects K.qnp
+            if self.backend.np.real(o) >= 0: # TODO: check for side effects K.qnp
                 r._eigenvalues = o * self._eigenvalues
             elif not self.backend.issparse(self.matrix):
                 r._eigenvalues = o * self._eigenvalues[::-1]
         if self._eigenvectors is not None:
-            if self.backend.cast(o).real > 0: # TODO: see above
+            if self.backend.np.real(o) > 0: # TODO: see above
                 r._eigenvectors = self._eigenvectors
             elif o == 0:
                 r._eigenvectors = self.eye(int(self._eigenvectors.shape[0]))
diff --git a/src/qibo/models/circuit.py b/src/qibo/models/circuit.py
index aa0045f91b..9897eaafa9 100644
--- a/src/qibo/models/circuit.py
+++ b/src/qibo/models/circuit.py
@@ -133,6 +133,7 @@ def __init__(self, nqubits, accelerators=None, density_matrix=False):
         self.measurement_gate_result = None
 
         self._final_state = None
+        self.compiled = None
         self.repeated_execution = False
 
         self.density_matrix = density_matrix
@@ -811,14 +812,37 @@ def final_state(self):
                                       "circuit is executed.")
         return self._final_state
 
+    def compile(self, backend=None):
+        if self.compiled:
+            raise_error(RuntimeError, "Circuit is already compiled.")
+        if not self.queue:
+            raise_error(RuntimeError, "Cannot compile circuit without gates.")
+        for gate in self.queue:
+            if isinstance(gate, gates.CallbackGate):
+                raise_error(NotImplementedError, "Circuit compilation is not available with callbacks.")
+        if backend is None:
+            from qibo.backends import GlobalBackend
+            backend = GlobalBackend()
+
+        from qibo.states import CircuitResult
+        executor = lambda state, nshots: backend.execute_circuit(self, state, nshots, return_array=True)
+        self.compiled = type("CompiledExecutor", (), {})()
+        self.compiled.executor = backend.compile(executor)
+        self.compiled.result = lambda state, nshots: CircuitResult(backend, self, state, nshots)
+
     def execute(self, initial_state=None, nshots=None):
         """Executes the circuit. Exact implementation depends on the backend.
 
         See :meth:`qibo.core.circuit.Circuit.execute` for more
         details.
         """
-        from qibo.backends import GlobalBackend
-        return GlobalBackend().execute_circuit(self, initial_state, nshots)
+        if self.compiled:
+            state = self.compiled.executor(initial_state, nshots)
+            self._final_state = self.compiled.result(state, nshots)
+            return self._final_state
+        else:
+            from qibo.backends import GlobalBackend
+            return GlobalBackend().execute_circuit(self, initial_state, nshots)
 
     def __call__(self, initial_state=None, nshots=None):
         """Equivalent to ``circuit.execute``."""
diff --git a/src/qibo/optimizers.py b/src/qibo/optimizers.py
index 1211c9617e..f17f6e8402 100644
--- a/src/qibo/optimizers.py
+++ b/src/qibo/optimizers.py
@@ -131,7 +131,7 @@ def newtonian(loss, initial_parameters, args=(), method='Powell',
     """
     if method == 'parallel_L-BFGS-B':  # pragma: no cover
         from qibo.parallel import _check_parallel_configuration
-        _check_parallel_configuration(processes)
+        _check_parallel_configuration(processes)  # pylint: disable=E1120
         o = ParallelBFGS(loss, args=args, processes=processes,
                          bounds=bounds, callback=callback, options=options)
         m = o.run(initial_parameters)
diff --git a/src/qibo/tests/conftest.py b/src/qibo/tests/conftest.py
index 539f46e682..bbf075d83f 100644
--- a/src/qibo/tests/conftest.py
+++ b/src/qibo/tests/conftest.py
@@ -12,7 +12,6 @@
     "qibo.tests.test_backends_agreement",
     "qibo.tests.test_backends_init",
     "qibo.tests.test_backends_matrices",
-    "qibo.tests.test_core_circuit_backpropagation",
     "qibo.tests.test_core_distcircuit_execution",
     "qibo.tests.test_core_distcircuit",
     "qibo.tests.test_core_distutils",
@@ -22,6 +21,7 @@
     "qibo.tests.test_models_evolution",
     "qibo.tests.test_models_qgan",
     "qibo.tests.test_models_variational",
+    "qibo.tests.test_parallel"
 }
 
 # backends to be tested
diff --git a/src/qibo/tests/test_callbacks.py b/src/qibo/tests/test_callbacks.py
index 4d5700eab3..cbe0b2af8e 100644
--- a/src/qibo/tests/test_callbacks.py
+++ b/src/qibo/tests/test_callbacks.py
@@ -161,23 +161,6 @@ def test_entropy_in_distributed_circuit(backend, accelerators, gateconf, target_
     backend.assert_allclose(values, target_entropy, atol=_atol)
 
 
-@pytest.mark.skip
-def test_entropy_in_compiled_circuit(backend):
-    """Check that entropy calculation works when circuit is compiled."""
-    from qibo import get_backend
-    entropy = callbacks.EntanglementEntropy([0])
-    c = Circuit(2)
-    c.add(gates.CallbackGate(entropy))
-    c.add(gates.H(0))
-    c.add(gates.CallbackGate(entropy))
-    c.add(gates.CNOT(0, 1))
-    c.add(gates.CallbackGate(entropy))
-    c.compile()
-    final_state = backend.execute_circuit(c)
-    values = [backend.to_numpy(x) for x in entropy[:]]
-    backend.assert_allclose(values, [0, 0, 1.0], atol=_atol)
-
-
 def test_entropy_multiple_executions(backend, accelerators):
     """Check entropy calculation when the callback is used in multiple executions."""
     target_c = Circuit(4)
@@ -300,7 +283,7 @@ def test_state_callback(backend, density_matrix, copy):
 
     target_state0 = np.array([1, 0, 1, 0]) / np.sqrt(2)
     target_state1 = np.ones(4) / 2.0
-    if not copy and str(backend) == "qibojit (numba)":
+    if not copy and backend.name == "qibojit":
         # when copy is disabled in the callback and in-place updates are used
         target_state0 = target_state1
     if density_matrix:
diff --git a/src/qibo/tests/test_core_circuit_backpropagation.py b/src/qibo/tests/test_core_circuit_backpropagation.py
deleted file mode 100644
index 479b60ec43..0000000000
--- a/src/qibo/tests/test_core_circuit_backpropagation.py
+++ /dev/null
@@ -1,50 +0,0 @@
-"""Tests Tensorflow's backpropagation when using `tf.Variable` parameters."""
-import numpy as np
-import pytest
-from qibo import K, gates
-from qibo.models import Circuit
-
-
-def test_variable_backpropagation(backend):
-    if not K.supports_gradients:
-        pytest.skip("Backpropagation is not supported by {}.".format(K.name))
-
-    theta = K.optimization.Variable(0.1234, dtype=K.dtypes('DTYPE'))
-    # TODO: Fix parametrized gates so that `Circuit` can be defined outside
-    # of the gradient tape
-    with K.optimization.GradientTape() as tape:
-        c = Circuit(1)
-        c.add(gates.X(0))
-        c.add(gates.RZ(0, theta))
-        loss = K.real(c()[-1])
-    grad = tape.gradient(loss, theta)
-
-    target_loss = np.cos(theta / 2.0)
-    K.assert_allclose(loss, target_loss)
-
-    target_grad = - np.sin(theta / 2.0) / 2.0
-    K.assert_allclose(grad, target_grad)
-
-
-def test_two_variables_backpropagation(backend):
-    if not K.supports_gradients:
-        pytest.skip("Backpropagation is not supported by {}.".format(K.name))
-
-    theta = K.optimization.Variable([0.1234, 0.4321], dtype=K.dtypes('DTYPE'))
-    # TODO: Fix parametrized gates so that `Circuit` can be defined outside
-    # of the gradient tape
-    with K.optimization.GradientTape() as tape:
-        c = Circuit(2)
-        c.add(gates.RX(0, theta[0]))
-        c.add(gates.RY(1, theta[1]))
-        loss = K.real(c()[0])
-    grad = tape.gradient(loss, theta)
-
-    t = np.array([0.1234, 0.4321]) / 2.0
-    target_loss = np.cos(t[0]) * np.cos(t[1])
-    K.assert_allclose(loss, target_loss)
-
-    target_grad1 = - np.sin(t[0]) * np.cos(t[1])
-    target_grad2 = - np.cos(t[0]) * np.sin(t[1])
-    target_grad = np.array([target_grad1, target_grad2]) / 2.0
-    K.assert_allclose(grad, target_grad)
diff --git a/src/qibo/tests/test_gates_density_matrix.py b/src/qibo/tests/test_gates_density_matrix.py
index 4e804c74b2..7795fb9646 100644
--- a/src/qibo/tests/test_gates_density_matrix.py
+++ b/src/qibo/tests/test_gates_density_matrix.py
@@ -64,7 +64,7 @@ def test_one_qubit_gates(backend, gatename, gatekwargs):
     gate = getattr(gates, gatename)(0, **gatekwargs)
     final_rho = apply_gates(backend, [gate], 1, initial_rho)
 
-    matrix = gate.asmatrix(backend)
+    matrix = backend.to_numpy(gate.asmatrix(backend))
     target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T)
     backend.assert_allclose(final_rho, target_rho)
 
@@ -75,7 +75,7 @@ def test_controlled_by_one_qubit_gates(backend, gatename):
     gate = getattr(gates, gatename)(1).controlled_by(0)
     final_rho = apply_gates(backend, [gate], 2, initial_rho)
 
-    matrix = backend.asmatrix(getattr(gates, gatename)(1))
+    matrix = backend.to_numpy(backend.asmatrix(getattr(gates, gatename)(1)))
     cmatrix = np.eye(4, dtype=matrix.dtype)
     cmatrix[2:, 2:] = matrix
     target_rho = np.einsum("ab,bc,cd->ad", cmatrix, initial_rho, cmatrix.conj().T)
@@ -95,7 +95,7 @@ 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.asmatrix(backend)
+    matrix = backend.to_numpy(gate.asmatrix(backend))
     target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T)
     backend.assert_allclose(final_rho, target_rho, atol=_atol)
 
@@ -106,7 +106,7 @@ def test_toffoli_gate(backend):
     gate = gates.TOFFOLI(0, 1, 2)
     final_rho = apply_gates(backend, [gate], 3, initial_rho)
 
-    matrix = gate.asmatrix(backend)
+    matrix = backend.to_numpy(gate.asmatrix(backend))
     target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T)
     backend.assert_allclose(final_rho, target_rho)
 
diff --git a/src/qibo/tests/test_hamiltonians.py b/src/qibo/tests/test_hamiltonians.py
index 92d38e65c6..e81a5f9454 100644
--- a/src/qibo/tests/test_hamiltonians.py
+++ b/src/qibo/tests/test_hamiltonians.py
@@ -15,6 +15,7 @@ def test_hamiltonian_init(backend):
     with pytest.raises(ValueError):
         H3 = hamiltonians.Hamiltonian(4, np.eye(10), backend=backend)
 
+
 @pytest.mark.parametrize("dtype", [np.int, np.float, np.complex, np.int32,
                                    np.int64, np.float32, np.float64,
                                    np.complex64, np.complex128])
@@ -47,6 +48,8 @@ def transformation_d(a, b, use_eye=False):
         H2 = hamiltonians.XXZ(nqubits=2, delta=1, backend=backend)
         mH1, mH2 = backend.to_numpy(H1.matrix), backend.to_numpy(H2.matrix)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         mH1 = random_sparse_matrix(backend, 64, sparse_type=sparse_type)
         mH2 = random_sparse_matrix(backend, 64, sparse_type=sparse_type)
         H1 = hamiltonians.Hamiltonian(6, mH1, backend=backend)
@@ -74,6 +77,8 @@ def test_hamiltonian_addition(backend, sparse_type):
         H1 = hamiltonians.Y(nqubits=3, backend=backend)
         H2 = hamiltonians.TFIM(nqubits=3, h=1.0, backend=backend)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         H1 = hamiltonians.Hamiltonian(6,random_sparse_matrix(backend, 64, sparse_type=sparse_type),
                                       backend=backend)
         H2 = hamiltonians.Hamiltonian(6, random_sparse_matrix(backend, 64, sparse_type=sparse_type),
@@ -117,6 +122,8 @@ def test_hamiltonian_matmul(backend, sparse_type):
         H1 = hamiltonians.TFIM(nqubits, h=1.0, backend=backend)
         H2 = hamiltonians.Y(nqubits, backend=backend)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         nqubits = 5
         nstates = 2 ** nqubits
         H1 = hamiltonians.Hamiltonian(nqubits, random_sparse_matrix(backend, nstates, sparse_type),
@@ -146,10 +153,12 @@ def test_hamiltonian_matmul_states(backend, sparse_type):
         nqubits = 3
         H = hamiltonians.TFIM(nqubits, h=1.0, backend=backend)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         nqubits = 5
         nstates = 2 ** nqubits
-        H = hamiltonians.Hamiltonian(nqubits, random_sparse_matrix(backend, nstates, sparse_type),
-                                     backend=backend)
+        matrix = random_sparse_matrix(backend, nstates, sparse_type)
+        H = hamiltonians.Hamiltonian(nqubits, matrix, backend=backend)
 
     hm = backend.to_numpy(H.matrix)
     v = random_complex(2 ** nqubits, dtype=hm.dtype)
@@ -173,6 +182,8 @@ def test_hamiltonian_expectation(backend, dense, density_matrix, sparse_type):
     if sparse_type is None:
         h = hamiltonians.XXZ(nqubits=3, delta=0.5, dense=dense, backend=backend)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         h = hamiltonians.Hamiltonian(6, random_sparse_matrix(backend, 64, sparse_type), backend=backend)
 
     matrix = backend.to_numpy(h.matrix)
@@ -209,6 +220,8 @@ def test_hamiltonian_eigenvalues(backend, dtype, sparse_type, dense):
     if sparse_type is None:
         H1 = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense, backend=backend)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         from scipy import sparse
         H1 = hamiltonians.XXZ(nqubits=5, delta=0.5, backend=backend)
         m = getattr(sparse, f"{sparse_type}_matrix")(backend.to_numpy(H1.matrix))
@@ -272,6 +285,8 @@ def test_hamiltonian_ground_state(backend, sparse_type, dense):
     if sparse_type is None:
         H = hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense, backend=backend)
     else:
+        if backend.name == "tensorflow":
+            pytest.skip("Tensorflow does not support operations with sparse matrices.")
         from scipy import sparse
         H = hamiltonians.XXZ(nqubits=5, delta=0.5, backend=backend)
         m = getattr(sparse, f"{sparse_type}_matrix")(backend.to_numpy(H.matrix))
@@ -291,6 +306,8 @@ def construct_hamiltonian():
         if sparse_type is None:
             return hamiltonians.XXZ(nqubits=2, delta=0.5, dense=dense, backend=backend)
         else:
+            if backend.name == "tensorflow":
+                pytest.skip("Tensorflow does not support operations with sparse matrices.")
             from scipy import sparse
             ham = hamiltonians.XXZ(nqubits=5, delta=0.5, backend=backend)
             m = getattr(sparse, f"{sparse_type}_matrix")(backend.to_numpy(ham.matrix))
diff --git a/src/qibo/tests/test_measurements.py b/src/qibo/tests/test_measurements.py
index 5e719e0e39..48cbb44be4 100644
--- a/src/qibo/tests/test_measurements.py
+++ b/src/qibo/tests/test_measurements.py
@@ -229,17 +229,13 @@ def test_circuit_copy_with_measurements(backend, accelerators):
         assert rg1[k] == rg2[k]
 
 
-@pytest.mark.skip
 def test_measurement_compiled_circuit(backend):
-    if backend.is_custom:
-        # use native gates because custom gates do not support compilation
-        pytest.skip("Custom backend does not support compilation.")
     c = models.Circuit(2)
     c.add(gates.X(0))
     c.add(gates.M(0))
     c.add(gates.M(1))
-    c.compile()
-    result = backend.execute_circuit(c, nshots=100)
+    c.compile(backend)
+    result = c(nshots=100)
     target_binary_samples = np.zeros((100, 2))
     target_binary_samples[:, 0] = 1
     assert_result(backend, result, 
diff --git a/src/qibo/tests/test_measurements_probabilistic.py b/src/qibo/tests/test_measurements_probabilistic.py
index 45333bd61d..2de8d4077e 100644
--- a/src/qibo/tests/test_measurements_probabilistic.py
+++ b/src/qibo/tests/test_measurements_probabilistic.py
@@ -60,20 +60,19 @@ def test_measurements_with_probabilistic_noise(backend):
     result = backend.execute_circuit(c, nshots=20)
     samples = result.samples()
 
-    np.random.seed(123)
     backend.set_seed(123)
     target_samples = []
     for _ in range(20):
         noiseless_c = models.Circuit(5)
         noiseless_c.add((gates.RX(i, t) for i, t in enumerate(thetas)))
         for i in range(5):
-            if np.random.random() < 0.2:
+            if backend.np.random.random() < 0.2:
                 noiseless_c.add(gates.Y(i))
-            if np.random.random() < 0.4:
+            if backend.np.random.random() < 0.4:
                 noiseless_c.add(gates.Z(i))
         noiseless_c.add(gates.M(*range(5)))
         result = backend.execute_circuit(noiseless_c, nshots=1)
-        target_samples.append(result.samples())
+        target_samples.append(backend.to_numpy(result.samples()))
     target_samples = np.concatenate(target_samples, axis=0)
     backend.assert_allclose(samples, target_samples)
 
diff --git a/src/qibo/tests/test_models_circuit_backpropagation.py b/src/qibo/tests/test_models_circuit_backpropagation.py
new file mode 100644
index 0000000000..182ec56c3b
--- /dev/null
+++ b/src/qibo/tests/test_models_circuit_backpropagation.py
@@ -0,0 +1,59 @@
+"""Tests Tensorflow's backpropagation when using `tf.Variable` parameters."""
+import pytest
+import numpy as np
+from qibo import gates
+from qibo.models import Circuit
+
+
+def construct_tensorflow_backend():
+    try:
+        from qibo.backends import construct_backend
+        backend = construct_backend("tensorflow")
+    except ModuleNotFoundError: # pragma: no cover
+        pytest.skip("Skipping backpropagation test because tensorflow is not available.")
+    return backend
+
+
+def test_variable_backpropagation():
+    backend = construct_tensorflow_backend()
+    import tensorflow as tf
+    theta = tf.Variable(0.1234, dtype="float64")
+    # TODO: Fix parametrized gates so that `Circuit` can be defined outside
+    # of the gradient tape
+    with tf.GradientTape() as tape:
+        c = Circuit(1)
+        c.add(gates.X(0))
+        c.add(gates.RZ(0, theta))
+        result = backend.execute_circuit(c)
+        loss = tf.math.real(result.state()[-1])
+    grad = tape.gradient(loss, theta)
+
+    target_loss = np.cos(theta / 2.0)
+    backend.assert_allclose(loss, target_loss)
+
+    target_grad = - np.sin(theta / 2.0) / 2.0
+    backend.assert_allclose(grad, target_grad)
+
+
+def test_two_variables_backpropagation():
+    backend = construct_tensorflow_backend()
+    import tensorflow as tf
+    theta = tf.Variable([0.1234, 0.4321], dtype="float64")
+    # TODO: Fix parametrized gates so that `Circuit` can be defined outside
+    # of the gradient tape
+    with tf.GradientTape() as tape:
+        c = Circuit(2)
+        c.add(gates.RX(0, theta[0]))
+        c.add(gates.RY(1, theta[1]))
+        result = backend.execute_circuit(c)
+        loss = tf.math.real(result.state()[0])
+    grad = tape.gradient(loss, theta)
+
+    t = np.array([0.1234, 0.4321]) / 2.0
+    target_loss = np.cos(t[0]) * np.cos(t[1])
+    backend.assert_allclose(loss, target_loss)
+
+    target_grad1 = - np.sin(t[0]) * np.cos(t[1])
+    target_grad2 = - np.cos(t[0]) * np.sin(t[1])
+    target_grad = np.array([target_grad1, target_grad2]) / 2.0
+    backend.assert_allclose(grad, target_grad)
diff --git a/src/qibo/tests/test_models_circuit_execution.py b/src/qibo/tests/test_models_circuit_execution.py
index 919547da96..19fb04ed3f 100644
--- a/src/qibo/tests/test_models_circuit_execution.py
+++ b/src/qibo/tests/test_models_circuit_execution.py
@@ -12,7 +12,6 @@ def test_eager_execute(backend, accelerators):
     backend.assert_allclose(final_state, target_state)
 
 
-@pytest.mark.skip
 def test_compiled_execute(backend):
     def create_circuit(theta = 0.1234):
         c = Circuit(2)
@@ -28,17 +27,15 @@ def create_circuit(theta = 0.1234):
 
     # Run eager circuit
     c1 = create_circuit()
-    r1 = c1.execute()
+    r1 = backend.execute_circuit(c1)
 
     # Run compiled circuit
     c2 = create_circuit()
-    c2.compile()
+    c2.compile(backend)
     r2 = c2()
-    init_state = c2.get_initial_state()
     np.testing.assert_allclose(r1, r2)
 
 
-@pytest.mark.skip
 def test_compiling_twice_exception(backend):
     """Check that compiling a circuit a second time raises error."""
     c = Circuit(2)
diff --git a/src/qibo/tests/test_models_circuit_features.py b/src/qibo/tests/test_models_circuit_features.py
index e5f906a1fd..beab25be10 100644
--- a/src/qibo/tests/test_models_circuit_features.py
+++ b/src/qibo/tests/test_models_circuit_features.py
@@ -298,11 +298,11 @@ def test_repeated_execute_pauli_noise_channel(backend):
         noiseless_c = Circuit(4)
         noiseless_c.add((gates.RY(i, t) for i, t in enumerate(thetas)))
         for i in range(4):
-            if np.random.random() < 0.1:
+            if backend.np.random.random() < 0.1:
                 noiseless_c.add(gates.X(i))
-            if np.random.random() < 0.2:
+            if backend.np.random.random() < 0.2:
                 noiseless_c.add(gates.Y(i))
-            if np.random.random() < 0.3:
+            if backend.np.random.random() < 0.3:
                 noiseless_c.add(gates.Z(i))
         result = backend.execute_circuit(noiseless_c)
         target_state.append(result.state(numpy=True))
@@ -319,15 +319,15 @@ def test_repeated_execute_with_noise(backend):
     backend.set_seed(1234)
     final_state = backend.execute_circuit(noisy_c, nshots=20)
 
-    np.random.seed(1234)
+    backend.set_seed(1234)
     target_state = []
     for _ in range(20):
         noiseless_c = Circuit(4)
         for i, t in enumerate(thetas):
             noiseless_c.add(gates.RY(i, theta=t))
-            if np.random.random() < 0.2:
+            if backend.np.random.random() < 0.2:
                 noiseless_c.add(gates.X(i))
-            if np.random.random() < 0.1:
+            if backend.np.random.random() < 0.1:
                 noiseless_c.add(gates.Z(i))
         result = backend.execute_circuit(noiseless_c)
         target_state.append(result.state(numpy=True))
diff --git a/src/qibo/tests/test_models_circuit_parametrized.py b/src/qibo/tests/test_models_circuit_parametrized.py
index 47f99a3ea9..52a470bd63 100644
--- a/src/qibo/tests/test_models_circuit_parametrized.py
+++ b/src/qibo/tests/test_models_circuit_parametrized.py
@@ -243,21 +243,24 @@ def test_set_parameters_with_double_variationallayer(backend, nqubits, trainable
     backend.assert_circuitclose(c, target_c)
 
 
-@pytest.mark.skip
-def test_variable_theta(backend):
+def test_variable_theta():
     """Check that parametrized gates accept `tf.Variable` parameters."""
-    backend = qibo.get_backend()
-    if backend != "tensorflow":
-        pytest.skip("Numpy backends do not support variable parameters.")
-    theta1 = K.optimization.Variable(0.1234, dtype=K.dtypes('DTYPE'))
-    theta2 = K.optimization.Variable(0.4321, dtype=K.dtypes('DTYPE'))
+    try:
+        from qibo.backends import construct_backend
+        backend = construct_backend("tensorflow")
+    except ModuleNotFoundError: # pragma: no cover
+        pytest.skip("Skipping variable test because tensorflow is not available.")
+    
+    import tensorflow as tf
+    theta1 = tf.Variable(0.1234, dtype="float64")
+    theta2 = tf.Variable(0.4321, dtype="float64")
     cvar = Circuit(2)
     cvar.add(gates.RX(0, theta1))
     cvar.add(gates.RY(1, theta2))
-    final_state = cvar()
+    final_state = backend.execute_circuit(cvar)
 
     c = Circuit(2)
     c.add(gates.RX(0, 0.1234))
     c.add(gates.RY(1, 0.4321))
-    target_state = c()
-    K.assert_allclose(final_state, target_state)
+    target_state = backend.execute_circuit(c)
+    backend.assert_allclose(final_state, target_state)
diff --git a/src/qibo/tests/utils.py b/src/qibo/tests/utils.py
index a010ed3534..85601d00f2 100644
--- a/src/qibo/tests/utils.py
+++ b/src/qibo/tests/utils.py
@@ -29,13 +29,14 @@ def random_density_matrix(nqubits):
     rho[ids, ids] = rho[ids, ids] / np.trace(rho)
     return rho
 
+
 def random_sparse_matrix(backend, n, sparse_type=None):
     if backend.name == "tensorflow":
         nonzero = int(0.1 * n * n)
         indices = np.random.randint(0, n, size=(nonzero, 2))
         data = np.random.random(nonzero) + 1j * np.random.random(nonzero)
         data = backend.cast(data)
-        return backend.sparse.SparseTensor(indices, data, (n, n))
+        return backend.tf.sparse.SparseTensor(indices, data, (n, n))
     else:
         re = sparse.rand(n, n, format=sparse_type)
         im = sparse.rand(n, n, format=sparse_type)