diff --git a/tests/test_backends_torch_gradients.py b/tests/test_backends_torch_gradients.py
index 2549e6bd45..7895c47ce8 100644
--- a/tests/test_backends_torch_gradients.py
+++ b/tests/test_backends_torch_gradients.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from qibo import gates, models
+from qibo import Circuit, gates
 from qibo.backends import PyTorchBackend, construct_backend
 from qibo.quantum_info import infidelity
 
@@ -16,21 +16,21 @@ def test_torch_gradients():
     target_state = backend.np.rand(4, dtype=backend.np.complex128)
     target_state = target_state / backend.np.norm(target_state)
     params = backend.np.rand(4, dtype=backend.np.float64, requires_grad=True)
-    c = models.Circuit(2)
-    c.add(gates.RX(0, params[0]))
-    c.add(gates.RY(1, params[1]))
-    c.add(gates.U2(1, params[2], params[3]))
+    circuit = Circuit(2)
+    circuit.add(gates.RX(0, params[0]))
+    circuit.add(gates.RY(1, params[1]))
+    circuit.add(gates.U2(1, params[2], params[3]))
 
     initial_params = params.clone()
     initial_loss = infidelity(
-        target_state, backend.execute_circuit(c).state(), backend=backend
+        target_state, backend.execute_circuit(circuit).state(), backend=backend
     )
 
     optimizer = optimizer([params], lr=0.01)
     for _ in range(nepochs):
         optimizer.zero_grad()
-        c.set_parameters(params)
-        final_state = backend.execute_circuit(c).state()
+        circuit.set_parameters(params)
+        final_state = backend.execute_circuit(circuit).state()
         loss = infidelity(target_state, final_state, backend=backend)
         loss.backward()
         grad = params.grad.clone().norm()
@@ -51,13 +51,13 @@ def test_torch_tensorflow_gradients():
 
     target_state = backend.np.tensor([0.0, 1.0], dtype=backend.np.complex128)
     param = backend.np.tensor([0.1], dtype=backend.np.float64, requires_grad=True)
-    c = models.Circuit(1)
-    c.add(gates.RX(0, param[0]))
+    circuit = Circuit(1)
+    circuit.add(gates.RX(0, param[0]))
 
     optimizer = backend.np.optim.SGD
     optimizer = optimizer([param], lr=1)
     c.set_parameters(param)
-    final_state = backend.execute_circuit(c).state()
+    final_state = backend.execute_circuit(circuit).state()
     loss = infidelity(target_state, final_state, backend=backend)
     loss.backward()
     torch_param_grad = param.grad.clone().item()
@@ -66,14 +66,14 @@ def test_torch_tensorflow_gradients():
 
     target_state = tf_backend.tf.constant([0.0, 1.0], dtype=tf_backend.tf.complex128)
     param = tf_backend.tf.Variable([0.1], dtype=tf_backend.tf.float64)
-    c = models.Circuit(1)
-    c.add(gates.RX(0, param[0]))
+    c = Circuit(1)
+    circuit.add(gates.RX(0, param[0]))
 
     optimizer = tf_backend.tf.optimizers.SGD(learning_rate=1.0)
 
     with tf_backend.tf.GradientTape() as tape:
         c.set_parameters(param)
-        final_state = tf_backend.execute_circuit(c).state()
+        final_state = tf_backend.execute_circuit(circuit).state()
         loss = infidelity(target_state, final_state, backend=tf_backend)
 
     grads = tape.gradient(loss, [param])
diff --git a/tests/test_hamiltonians_terms.py b/tests/test_hamiltonians_terms.py
index 6d3dafdbad..0addb08adc 100644
--- a/tests/test_hamiltonians_terms.py
+++ b/tests/test_hamiltonians_terms.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from qibo import gates, matrices, models
+from qibo import Circuit, gates, matrices
 from qibo.hamiltonians import terms
 from qibo.quantum_info import random_density_matrix, random_statevector
 
@@ -49,15 +49,17 @@ def test_hamiltonian_term_gates(backend):
 
     initial_state = random_statevector(2**nqubits, backend=backend)
     final_state = term(backend, backend.np.copy(initial_state), nqubits)
-    c = models.Circuit(nqubits)
-    c.add(gates.Unitary(matrix, 1, 2))
-    target_state = backend.execute_circuit(c, backend.np.copy(initial_state)).state()
+    circuit = Circuit(nqubits)
+    circuit.add(gates.Unitary(matrix, 1, 2))
+    target_state = backend.execute_circuit(
+        circuit, backend.np.copy(initial_state)
+    ).state()
     backend.assert_allclose(final_state, target_state)
 
 
 def test_hamiltonian_term_exponentiation(backend):
     """Test exp gate application of ``HamiltonianTerm``."""
-    from scipy.linalg import expm
+    from scipy.linalg import expm  # pylint: disable=C0415
 
     matrix = np.random.random((2, 2))
     term = terms.HamiltonianTerm(matrix, 1)
diff --git a/tests/test_measurements_collapse.py b/tests/test_measurements_collapse.py
index 2944a67783..ce38b039a9 100644
--- a/tests/test_measurements_collapse.py
+++ b/tests/test_measurements_collapse.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from qibo import gates, models
+from qibo import Circuit, gates
 from qibo.quantum_info import random_density_matrix, random_statevector
 
 
@@ -13,7 +13,7 @@
 )
 def test_measurement_collapse(backend, nqubits, targets):
     initial_state = random_statevector(2**nqubits, backend=backend)
-    c = models.Circuit(nqubits)
+    c = Circuit(nqubits)
     for q in np.random.randint(nqubits, size=np.random.randint(nqubits, size=1)):
         c.add(gates.H(q))
     r = c.add(gates.M(*targets, collapse=True))
@@ -39,7 +39,7 @@ def assign_value(rho, index, value):
         return rho
 
     initial_rho = random_density_matrix(2**nqubits, backend=backend)
-    c = models.Circuit(nqubits, density_matrix=True)
+    c = Circuit(nqubits, density_matrix=True)
     r = c.add(gates.M(*targets, collapse=True))
     final_rho = backend.execute_circuit(c, backend.np.copy(initial_rho), nshots=1)
 
@@ -60,7 +60,7 @@ def assign_value(rho, index, value):
 
 
 def test_measurement_collapse_bitflip_noise(backend):
-    c = models.Circuit(4)
+    c = Circuit(4)
     with pytest.raises(NotImplementedError):
         output = c.add(gates.M(0, 1, p0=0.2, collapse=True))
 
@@ -68,7 +68,7 @@ def test_measurement_collapse_bitflip_noise(backend):
 @pytest.mark.parametrize("density_matrix", [True, False])
 @pytest.mark.parametrize("effect", [False, True])
 def test_measurement_result_parameters(backend, effect, density_matrix):
-    c = models.Circuit(4, density_matrix=density_matrix)
+    c = Circuit(4, density_matrix=density_matrix)
     if effect:
         c.add(gates.X(0))
     r = c.add(gates.M(0, collapse=True))
@@ -76,7 +76,7 @@ def test_measurement_result_parameters(backend, effect, density_matrix):
     if not density_matrix:
         c.add(gates.M(0))
 
-    target_c = models.Circuit(4, density_matrix=density_matrix)
+    target_c = Circuit(4, density_matrix=density_matrix)
     if effect:
         target_c.add(gates.X(0))
         target_c.add(gates.RX(1, theta=np.pi / 4))
@@ -94,7 +94,7 @@ def test_measurement_result_parameters(backend, effect, density_matrix):
 def test_measurement_result_parameters_random(backend):
     initial_state = random_density_matrix(2**4, backend=backend)
     backend.set_seed(123)
-    c = models.Circuit(4, density_matrix=True)
+    c = Circuit(4, density_matrix=True)
     r = c.add(gates.M(1, collapse=True))
     c.add(gates.RY(0, theta=np.pi * r.symbols[0] / 5))
     c.add(gates.RX(2, theta=np.pi * r.symbols[0] / 4))
@@ -103,13 +103,13 @@ def test_measurement_result_parameters_random(backend):
     )
 
     backend.set_seed(123)
-    c = models.Circuit(4, density_matrix=True)
+    c = Circuit(4, density_matrix=True)
     m = c.add(gates.M(1, collapse=True))
     target_state = backend.execute_circuit(
         c, initial_state=backend.np.copy(initial_state), nshots=1
     ).state()
     if int(m.symbols[0].outcome()):
-        c = models.Circuit(4, density_matrix=True)
+        c = Circuit(4, density_matrix=True)
         c.add(gates.RY(0, theta=np.pi / 5))
         c.add(gates.RX(2, theta=np.pi / 4))
         target_state = backend.execute_circuit(c, initial_state=target_state)
@@ -120,7 +120,7 @@ def test_measurement_result_parameters_random(backend):
 def test_measurement_result_parameters_repeated_execution(backend, use_loop):
     initial_state = random_density_matrix(2**4, backend=backend)
     backend.set_seed(123)
-    c = models.Circuit(4, density_matrix=True)
+    c = Circuit(4, density_matrix=True)
     r = c.add(gates.M(1, collapse=True))
     c.add(gates.RX(2, theta=np.pi * r.symbols[0] / 4))
     if use_loop:
@@ -139,7 +139,7 @@ def test_measurement_result_parameters_repeated_execution(backend, use_loop):
     backend.set_seed(123)
     target_states = []
     for _ in range(20):
-        c = models.Circuit(4, density_matrix=True)
+        c = Circuit(4, density_matrix=True)
         m = c.add(gates.M(1, collapse=True))
         target_state = backend.execute_circuit(
             c, backend.np.copy(initial_state), nshots=1
@@ -157,7 +157,7 @@ def test_measurement_result_parameters_repeated_execution(backend, use_loop):
 def test_measurement_result_parameters_repeated_execution_final_measurements(backend):
     initial_state = random_density_matrix(2**4, backend=backend)
     backend.set_seed(123)
-    c = models.Circuit(4, density_matrix=True)
+    c = Circuit(4, density_matrix=True)
     r = c.add(gates.M(1, collapse=True))
     c.add(gates.RY(0, theta=np.pi * r.symbols[0] / 3))
     c.add(gates.RY(2, theta=np.pi * r.symbols[0] / 4))
@@ -170,12 +170,12 @@ def test_measurement_result_parameters_repeated_execution_final_measurements(bac
     backend.set_seed(123)
     target_samples = []
     for _ in range(30):
-        c = models.Circuit(4, density_matrix=True)
+        c = Circuit(4, density_matrix=True)
         m = c.add(gates.M(1, collapse=True))
         target_state = backend.execute_circuit(
             c, backend.np.copy(initial_state), nshots=1
         ).state()
-        c = models.Circuit(4, density_matrix=True)
+        c = Circuit(4, density_matrix=True)
         if int(m.symbols[0].outcome()):
             c.add(gates.RY(0, theta=np.pi / 3))
             c.add(gates.RY(2, theta=np.pi / 4))
@@ -188,14 +188,14 @@ def test_measurement_result_parameters_repeated_execution_final_measurements(bac
 def test_measurement_result_parameters_multiple_qubits(backend):
     initial_state = random_density_matrix(2**4, backend=backend)
     backend.set_seed(123)
-    c = models.Circuit(4, density_matrix=True)
+    c = Circuit(4, density_matrix=True)
     r = c.add(gates.M(0, 1, 2, collapse=True))
     c.add(gates.RY(1, theta=np.pi * r.symbols[0] / 5))
     c.add(gates.RX(3, theta=np.pi * r.symbols[2] / 3))
     final_state = backend.execute_circuit(c, backend.np.copy(initial_state), nshots=1)
 
     backend.set_seed(123)
-    c = models.Circuit(4, density_matrix=True)
+    c = Circuit(4, density_matrix=True)
     m = c.add(gates.M(0, 1, 2, collapse=True))
     target_state = backend.execute_circuit(
         c, backend.np.copy(initial_state), nshots=1
@@ -217,7 +217,7 @@ def test_measurement_result_parameters_multiple_qubits(backend):
 @pytest.mark.parametrize("nqubits,targets", [(5, [2, 4]), (6, [3, 5])])
 def test_measurement_collapse_distributed(backend, accelerators, nqubits, targets):
     initial_state = random_density_matrix(2**nqubits, backend=backend)
-    c = models.Circuit(nqubits, accelerators, density_matrix=True)
+    c = Circuit(nqubits, accelerators, density_matrix=True)
     m = c.add(gates.M(*targets, collapse=True))
     result = backend.execute_circuit(c, np.copy(initial_state), nshots=1).state()
     slicer = 2 * nqubits * [slice(None)]
@@ -235,13 +235,13 @@ def test_measurement_collapse_distributed(backend, accelerators, nqubits, target
 
 def test_collapse_after_measurement(backend):
     qubits = [0, 2, 3]
-    c = models.Circuit(5, density_matrix=True)
+    c = Circuit(5, density_matrix=True)
     c.add(gates.H(i) for i in range(5))
     m = c.add(gates.M(*qubits, collapse=True))
     c.add(gates.H(i) for i in range(5))
     final_state = backend.execute_circuit(c, nshots=1)
 
-    ct = models.Circuit(5, density_matrix=True)
+    ct = Circuit(5, density_matrix=True)
     bitstring = [r.outcome() for r in m.symbols]
     for i, r in zip(qubits, bitstring):
         if r:
@@ -252,7 +252,7 @@ def test_collapse_after_measurement(backend):
 
 
 def test_collapse_error(backend):
-    c = models.Circuit(1)
+    c = Circuit(1)
     m = c.add(gates.M(0, collapse=True))
     with pytest.raises(Exception) as exc_info:
         backend.execute_circuit(c)
diff --git a/tests/test_measurements_probabilistic.py b/tests/test_measurements_probabilistic.py
index 742c4838b2..a831e23ba8 100644
--- a/tests/test_measurements_probabilistic.py
+++ b/tests/test_measurements_probabilistic.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from qibo import gates, models
+from qibo import Circuit, gates
 
 from .test_measurements import assert_result
 
@@ -12,11 +12,11 @@
 def test_probabilistic_measurement(backend, accelerators, use_samples):
     # set single-thread to fix the random values generated from the frequency custom op
     backend.set_threads(1)
-    c = models.Circuit(4, accelerators)
-    c.add(gates.H(0))
-    c.add(gates.H(1))
-    c.add(gates.M(0, 1))
-    result = backend.execute_circuit(c, nshots=1000)
+    circuit = Circuit(4, accelerators)
+    circuit.add(gates.H(0))
+    circuit.add(gates.H(1))
+    circuit.add(gates.M(0, 1))
+    result = backend.execute_circuit(circuit, nshots=1000)
 
     backend.set_seed(1234)
     if use_samples:
@@ -33,11 +33,11 @@ def test_probabilistic_measurement(backend, accelerators, use_samples):
 def test_sample_frequency_agreement(backend):
     # set single-thread to fix the random values generated from the frequency custom op
     backend.set_threads(1)
-    c = models.Circuit(2)
-    c.add(gates.H(0))
-    c.add(gates.H(1))
-    c.add(gates.M(0, 1))
-    result = backend.execute_circuit(c, nshots=1000)
+    circuit = Circuit(2)
+    circuit.add(gates.H(0))
+    circuit.add(gates.H(1))
+    circuit.add(gates.M(0, 1))
+    result = backend.execute_circuit(circuit, nshots=1000)
 
     backend.set_seed(1234)
     target_frequencies = result.frequencies(binary=False)
@@ -54,9 +54,9 @@ def test_unbalanced_probabilistic_measurement(backend, use_samples):
     # set single-thread to fix the random values generated from the frequency custom op
     backend.set_threads(1)
     state = np.array([1, 1, 1, np.sqrt(3)]) / np.sqrt(6)
-    c = models.Circuit(2)
-    c.add(gates.M(0, 1))
-    result = backend.execute_circuit(c, initial_state=np.copy(state), nshots=1000)
+    circuit = Circuit(2)
+    circuit.add(gates.M(0, 1))
+    result = backend.execute_circuit(circuit, initial_state=np.copy(state), nshots=1000)
 
     backend.set_seed(1234)
     if use_samples:
@@ -75,14 +75,14 @@ def test_unbalanced_probabilistic_measurement(backend, use_samples):
 def test_measurements_with_probabilistic_noise(backend):
     """Check measurements when simulating noise with repeated execution."""
     thetas = np.random.random(5)
-    c = models.Circuit(5)
-    c.add((gates.RX(i, t) for i, t in enumerate(thetas)))
-    c.add(
+    circuit = Circuit(5)
+    circuit.add((gates.RX(i, t) for i, t in enumerate(thetas)))
+    circuit.add(
         gates.PauliNoiseChannel(i, list(zip(["Y", "Z"], [0.2, 0.4]))) for i in range(5)
     )
-    c.add(gates.M(*range(5)))
+    circuit.add(gates.M(*range(5)))
     backend.set_seed(123)
-    result = backend.execute_circuit(c, nshots=20)
+    result = backend.execute_circuit(circuit, nshots=20)
     samples = result.samples()
 
     backend.set_seed(123)
@@ -90,14 +90,14 @@ def test_measurements_with_probabilistic_noise(backend):
     channel_gates = [gates.Y, gates.Z]
     probs = [0.2, 0.4, 0.4]
     for _ in range(20):
-        noiseless_c = models.Circuit(5)
-        noiseless_c.add((gates.RX(i, t) for i, t in enumerate(thetas)))
+        noiseless_circuit = Circuit(5)
+        noiseless_circuit.add((gates.RX(i, t) for i, t in enumerate(thetas)))
         for i in range(5):
             index = backend.sample_shots(probs, 1)[0]
             if index != len(channel_gates):
-                noiseless_c.add(channel_gates[index](i))
-        noiseless_c.add(gates.M(*range(5)))
-        result = backend.execute_circuit(noiseless_c, nshots=1)
+                noiseless_circuit.add(channel_gates[index](i))
+        noiseless_circuit.add(gates.M(*range(5)))
+        result = backend.execute_circuit(noiseless_circuit, nshots=1)
         target_samples.append(backend.to_numpy(result.samples()))
     target_samples = np.concatenate(target_samples, axis=0)
     backend.assert_allclose(samples, target_samples)
@@ -109,11 +109,11 @@ def test_measurements_with_probabilistic_noise(backend):
 def test_post_measurement_bitflips_on_circuit(backend, accelerators, i, probs):
     """Check bitflip errors on circuit measurements."""
     backend.set_seed(123)
-    c = models.Circuit(5, accelerators=accelerators)
-    c.add([gates.X(0), gates.X(2), gates.X(3)])
-    c.add(gates.M(0, 1, p0={0: probs[0], 1: probs[1]}))
-    c.add(gates.M(3, p0=probs[2]))
-    result = backend.execute_circuit(c, nshots=30)
+    circuit = Circuit(5, accelerators=accelerators)
+    circuit.add([gates.X(0), gates.X(2), gates.X(3)])
+    circuit.add(gates.M(0, 1, p0={0: probs[0], 1: probs[1]}))
+    circuit.add(gates.M(3, p0=probs[2]))
+    result = backend.execute_circuit(circuit, nshots=30)
     freqs = result.frequencies(binary=False)
     targets = backend._test_regressions("test_post_measurement_bitflips_on_circuit")
     assert freqs == targets[i]
@@ -123,11 +123,11 @@ def test_post_measurement_bitflips_on_circuit_result(backend):
     """Check bitflip errors on ``CircuitResult`` objects."""
     thetas = np.random.random(4)
     backend.set_seed(123)
-    c = models.Circuit(4)
-    c.add((gates.RX(i, theta=t) for i, t in enumerate(thetas)))
-    c.add(gates.M(0, 1, register_name="a", p0={0: 0.2, 1: 0.4}))
-    c.add(gates.M(3, register_name="b", p0=0.3))
-    result = backend.execute_circuit(c, nshots=30)
+    circuit = Circuit(4)
+    circuit.add((gates.RX(i, theta=t) for i, t in enumerate(thetas)))
+    circuit.add(gates.M(0, 1, register_name="a", p0={0: 0.2, 1: 0.4}))
+    circuit.add(gates.M(3, register_name="b", p0=0.3))
+    result = backend.execute_circuit(circuit, nshots=30)
     samples = result.samples(binary=True)
     register_samples = result.samples(binary=True, registers=True)
     backend.assert_allclose(register_samples["a"], samples[:, :2])
@@ -146,10 +146,10 @@ def test_post_measurement_bitflips_on_circuit_result(backend):
 def test_measurementresult_apply_bitflips(backend, i, p0, p1):
     from qibo.result import CircuitResult
 
-    c = models.Circuit(3)
-    c.add(gates.M(*range(3)))
+    circuit = Circuit(3)
+    circuit.add(gates.M(*range(3)))
     state = backend.zero_state(8)
-    result = CircuitResult(state, c.measurements, backend)
+    result = CircuitResult(state, circuit.measurements, backend)
     result._samples = backend.cast(np.zeros((10, 3)), dtype="int32")
     backend.set_seed(123)
     noisy_samples = result.apply_bitflips(p0, p1)
diff --git a/tests/test_models_variational.py b/tests/test_models_variational.py
index 8704aaa31f..1c49f4bacf 100644
--- a/tests/test_models_variational.py
+++ b/tests/test_models_variational.py
@@ -8,9 +8,11 @@
 import pytest
 from scipy.linalg import expm
 
-from qibo import gates, hamiltonians, models
+from qibo import Circuit, gates
+from qibo.hamiltonians import TFIM, XXZ, X, Y
 from qibo.models.utils import cvar, gibbs
-from qibo.quantum_info import random_statevector
+from qibo.models.variational import AAVQE, FALQON, QAOA, VQE
+from qibo.quantum_info.random_ensembles import random_statevector
 
 REGRESSION_FOLDER = pathlib.Path(__file__).with_name("regressions")
 
@@ -63,14 +65,14 @@ def myloss(parameters, circuit, target):
     nlayers = 4
 
     # Create variational circuit
-    c = models.Circuit(nqubits)
+    circuit = Circuit(nqubits)
     for _ in range(nlayers):
-        c.add(gates.RY(q, theta=0) for q in range(nqubits))
-        c.add(gates.CZ(q, q + 1) for q in range(0, nqubits - 1, 2))
-        c.add(gates.RY(q, theta=0) for q in range(nqubits))
-        c.add(gates.CZ(q, q + 1) for q in range(1, nqubits - 2, 2))
-        c.add(gates.CZ(0, nqubits - 1))
-    c.add(gates.RY(q, theta=0) for q in range(nqubits))
+        circuit.add(gates.RY(qubit, theta=0.0) for qubit in range(nqubits))
+        circuit.add(gates.CZ(qubit, qubit + 1) for qubit in range(0, nqubits - 1, 2))
+        circuit.add(gates.RY(qubit, theta=0.0) for qubit in range(nqubits))
+        circuit.add(gates.CZ(qubit, qubit + 1) for qubit in range(1, nqubits - 2, 2))
+        circuit.add(gates.CZ(0, nqubits - 1))
+    circuit.add(gates.RY(qubit, theta=0) for qubit in range(nqubits))
 
     # Optimize starting from a random guess for the variational parameters
     np.random.seed(0)
@@ -78,8 +80,13 @@ def myloss(parameters, circuit, target):
     data = np.random.normal(0, 1, size=2**nqubits)
 
     # perform optimization
-    best, params, _ = optimize(
-        myloss, x0, args=(c, data), method=method, options=options, compile=compile
+    _, params, _ = optimize(
+        myloss,
+        x0,
+        args=(circuit, data),
+        method=method,
+        options=options,
+        compile=compile,
     )
     if filename is not None:
         assert_regression_fixture(backend, params, filename)
@@ -113,7 +120,7 @@ def test_vqe(backend, method, options, compile, filename):
     backend.set_threads(1)
     nqubits = 3
     layers = 4
-    circuit = models.Circuit(nqubits)
+    circuit = Circuit(nqubits)
     for l in range(layers):
         for q in range(nqubits):
             circuit.add(gates.RY(q, theta=1.0))
@@ -126,14 +133,14 @@ def test_vqe(backend, method, options, compile, filename):
         circuit.add(gates.CZ(0, nqubits - 1))
     for q in range(nqubits):
         circuit.add(gates.RY(q, theta=1.0))
-    hamiltonian = hamiltonians.XXZ(nqubits=nqubits, backend=backend)
+    hamiltonian = XXZ(nqubits=nqubits, backend=backend)
     np.random.seed(0)
     initial_parameters = backend.cast(
         np.random.uniform(0, 2 * np.pi, 2 * nqubits * layers + nqubits), dtype="float64"
     )
     if backend.name == "pytorch":
         initial_parameters.requires_grad = True
-    v = models.VQE(circuit, hamiltonian)
+    v = VQE(circuit, hamiltonian)
 
     loss_values = []
 
@@ -176,8 +183,8 @@ def callback(parameters, loss_values=loss_values, vqe=v):
     ],
 )
 def test_qaoa_execution(backend, solver, dense, accel=None):
-    h = hamiltonians.TFIM(6, h=1.0, dense=dense, backend=backend)
-    m = hamiltonians.X(6, dense=dense, backend=backend)
+    h = TFIM(6, h=1.0, dense=dense, backend=backend)
+    m = X(6, dense=dense, backend=backend)
     # Trotter and RK require small p's!
     params = 0.01 * (1 - 2 * np.random.random(4))
     state = random_statevector(2**6, backend=backend)
@@ -199,7 +206,7 @@ def test_qaoa_execution(backend, solver, dense, accel=None):
             u = expm(-1j * p * h_matrix)
         target_state = backend.cast(u) @ target_state
 
-    qaoa = models.QAOA(h, mixer=m, solver=solver, accelerators=accel)
+    qaoa = QAOA(h, mixer=m, solver=solver, accelerators=accel)
     qaoa.set_parameters(params)
     final_state = qaoa(backend.cast(state, copy=True))
     backend.assert_allclose(final_state, target_state, atol=atol)
@@ -214,13 +221,13 @@ def test_qaoa_callbacks(backend, accelerators):
 
     # use ``Y`` Hamiltonian so that there are no errors
     # in the Trotter decomposition
-    h = hamiltonians.Y(5, backend=backend)
+    h = Y(5, backend=backend)
     energy = callbacks.Energy(h)
     params = 0.1 * np.random.random(4)
     state = random_statevector(2**5, backend=backend)
 
-    ham = hamiltonians.Y(5, dense=False, backend=backend)
-    qaoa = models.QAOA(ham, callbacks=[energy], accelerators=accelerators)
+    ham = Y(5, dense=False, backend=backend)
+    qaoa = QAOA(ham, callbacks=[energy], accelerators=accelerators)
     qaoa.set_parameters(params)
     final_state = qaoa(backend.cast(state, copy=True))
 
@@ -243,22 +250,22 @@ def test_qaoa_callbacks(backend, accelerators):
 def test_qaoa_errors(backend):
     # Invalid Hamiltonian type
     with pytest.raises(TypeError):
-        qaoa = models.QAOA("test")
+        qaoa = QAOA("test")
     # Hamiltonians of different type
-    h = hamiltonians.TFIM(4, h=1.0, dense=False, backend=backend)
-    m = hamiltonians.X(4, dense=True, backend=backend)
+    h = TFIM(4, h=1.0, dense=False, backend=backend)
+    m = X(4, dense=True, backend=backend)
     with pytest.raises(TypeError):
-        qaoa = models.QAOA(h, mixer=m)
+        qaoa = QAOA(h, mixer=m)
     # Hamiltonians acting on different qubit numbers
-    h = hamiltonians.TFIM(6, h=1.0, backend=backend)
-    m = hamiltonians.X(4, backend=backend)
+    h = TFIM(6, h=1.0, backend=backend)
+    m = X(4, backend=backend)
     with pytest.raises(ValueError):
-        qaoa = models.QAOA(h, mixer=m)
+        qaoa = QAOA(h, mixer=m)
     # distributed execution with RK solver
     with pytest.raises(NotImplementedError):
-        qaoa = models.QAOA(h, solver="rk4", accelerators={"/GPU:0": 2})
+        qaoa = QAOA(h, solver="rk4", accelerators={"/GPU:0": 2})
     # minimize with odd number of parameters
-    qaoa = models.QAOA(h)
+    qaoa = QAOA(h)
     with pytest.raises(ValueError):
         qaoa.minimize(np.random.random(5))
 
@@ -278,8 +285,8 @@ def test_qaoa_optimization(backend, method, options, dense, filename):
         pytest.skip("Skipping SGD test for unsupported backend.")
     if method != "sgd" and backend.name in ("tensorflow", "pytorch"):
         pytest.skip("Skipping scipy optimizers for tensorflow and pytorch.")
-    h = hamiltonians.XXZ(3, dense=dense, backend=backend)
-    qaoa = models.QAOA(h)
+    h = XXZ(3, dense=dense, backend=backend)
+    qaoa = QAOA(h)
     initial_p = backend.cast([0.05, 0.06, 0.07, 0.08], dtype="float64")
     if backend.name == "pytorch":
         initial_p.requires_grad = True
@@ -299,8 +306,8 @@ def test_qaoa_optimization(backend, method, options, dense, filename):
 
 @pytest.mark.parametrize(test_names, test_values)
 def test_falqon_optimization(backend, delta_t, max_layers, tolerance, filename):
-    h = hamiltonians.XXZ(3, backend=backend)
-    falqon = models.FALQON(h)
+    h = XXZ(3, backend=backend)
+    falqon = FALQON(h)
     best, params, extra = falqon.minimize(delta_t, max_layers, tol=tolerance)
     if filename is not None:
         assert_regression_fixture(backend, params, filename)
@@ -312,8 +319,8 @@ def __call__(self, x):
             return np.sum(x)
 
     callback = TestCallback()
-    h = hamiltonians.XXZ(3, backend=backend)
-    falqon = models.FALQON(h)
+    h = XXZ(3, backend=backend)
+    falqon = FALQON(h)
     best, params, extra = falqon.minimize(0.1, 5, callback=callback)
     assert len(extra["callbacks"]) == 5
 
@@ -332,7 +339,7 @@ def test_aavqe(backend, method, options, compile, filename):
 
     nqubits = 4
     layers = 1
-    circuit = models.Circuit(nqubits)
+    circuit = Circuit(nqubits)
 
     for l in range(layers):
         for q in range(nqubits):
@@ -347,12 +354,10 @@ def test_aavqe(backend, method, options, compile, filename):
     for q in range(nqubits):
         circuit.add(gates.RY(q, theta=1.0))
 
-    easy_hamiltonian = hamiltonians.X(nqubits, backend=backend)
-    problem_hamiltonian = hamiltonians.XXZ(nqubits, backend=backend)
+    easy_hamiltonian = X(nqubits, backend=backend)
+    problem_hamiltonian = XXZ(nqubits, backend=backend)
     s = lambda t: t
-    aavqe = models.AAVQE(
-        circuit, easy_hamiltonian, problem_hamiltonian, s, nsteps=10, t_max=1
-    )
+    aavqe = AAVQE(circuit, easy_hamiltonian, problem_hamiltonian, s, nsteps=10, t_max=1)
     np.random.seed(0)
     initial_parameters = np.random.uniform(0, 2 * np.pi, 2 * nqubits * layers + nqubits)
     best, params = aavqe.minimize(
@@ -374,8 +379,8 @@ def test_aavqe(backend, method, options, compile, filename):
 def test_custom_loss(test_input, test_param, expected):
     from qibo import hamiltonians
 
-    h = hamiltonians.XXZ(3)
-    qaoa = models.QAOA(h)
+    h = XXZ(3)
+    qaoa = QAOA(h)
     initial_p = [0.314, 0.22, 0.05, 0.59]
     best, params, _ = qaoa.minimize(
         initial_p, loss_func=test_input, loss_func_param=test_param
diff --git a/tests/test_result.py b/tests/test_result.py
index 7b37466a2f..3d57f0140a 100644
--- a/tests/test_result.py
+++ b/tests/test_result.py
@@ -3,26 +3,25 @@
 import numpy as np
 import pytest
 
-from qibo import Circuit, gates, models
-from qibo.config import raise_error
+from qibo import Circuit, gates
 from qibo.result import CircuitResult, MeasurementOutcomes, load_result
 
 
 @pytest.mark.parametrize("qubits", [None, [0], [1, 2]])
 def test_measurementoutcomes_probabilities(backend, qubits):
-    c = Circuit(3)
-    c.add(gates.H(0))
-    c.add(gates.M(0, 2))
-    global_probs = backend.execute_circuit(c).probabilities(qubits=[0, 2])
+    circuit = Circuit(3)
+    circuit.add(gates.H(0))
+    circuit.add(gates.M(0, 2))
+    global_probs = backend.execute_circuit(circuit).probabilities(qubits=[0, 2])
     probabilities = (
-        backend.execute_circuit(c).probabilities(qubits=qubits)
+        backend.execute_circuit(circuit).probabilities(qubits=qubits)
         if qubits is not None
-        else backend.execute_circuit(c).probabilities(qubits=[0, 2])
+        else backend.execute_circuit(circuit).probabilities(qubits=[0, 2])
     )
-    c.has_collapse = True
+    circuit.has_collapse = True
     if qubits is not None and 1 in qubits:
         with pytest.raises(RuntimeError) as excinfo:
-            repeated_probabilities = backend.execute_circuit(c).probabilities(
+            repeated_probabilities = backend.execute_circuit(circuit).probabilities(
                 qubits=qubits
             )
             assert (
@@ -30,28 +29,28 @@ def test_measurementoutcomes_probabilities(backend, qubits):
                 == f"Asking probabilities for qubits {qubits}, but only qubits [0,2] were measured."
             )
     else:
-        repeated_probabilities = backend.execute_circuit(c, nshots=1000).probabilities(
-            qubits=qubits
-        )
+        repeated_probabilities = backend.execute_circuit(
+            circuit, nshots=1000
+        ).probabilities(qubits=qubits)
         result = MeasurementOutcomes(
-            c.measurements, backend=backend, probabilities=global_probs
+            circuit.measurements, backend=backend, probabilities=global_probs
         )
         backend.assert_allclose(probabilities, repeated_probabilities, atol=1e-1)
         backend.assert_allclose(result.probabilities(qubits), probabilities, atol=1e-1)
 
 
 def test_measurementoutcomes_samples_from_measurements(backend):
-    c = Circuit(3)
-    c.add(gates.H(0))
-    c.add(gates.M(0, 2))
-    res = backend.execute_circuit(c)
+    circuit = Circuit(3)
+    circuit.add(gates.H(0))
+    circuit.add(gates.M(0, 2))
+    res = backend.execute_circuit(circuit)
     samples = res.samples()
-    outcome = MeasurementOutcomes(c.measurements, backend=backend)
+    outcome = MeasurementOutcomes(circuit.measurements, backend=backend)
     backend.assert_allclose(samples, outcome.samples())
 
 
 def test_circuit_result_error(backend):
-    c = models.Circuit(1)
+    c = Circuit(1)
     state = np.array([1, 0])
     with pytest.raises(Exception) as exc_info:
         CircuitResult(state, c.measurements, backend)
@@ -62,7 +61,7 @@ def test_circuit_result_error(backend):
 
 
 def test_measurement_gate_dump_load(backend):
-    c = models.Circuit(2)
+    c = Circuit(2)
     c.add(gates.M(1, 0, basis=[gates.Z, gates.X]))
     m = c.measurements
     load = m[0].to_json()
@@ -72,7 +71,7 @@ def test_measurement_gate_dump_load(backend):
 
 @pytest.mark.parametrize("agnostic_load", [False, True])
 def test_measurementoutcomes_dump_load(backend, agnostic_load):
-    c = models.Circuit(2)
+    c = Circuit(2)
     c.add(gates.M(1, 0, basis=[gates.Z, gates.X]))
     # just to trigger repeated execution and test MeasurementOutcomes
     c.has_collapse = True
@@ -91,7 +90,7 @@ def test_measurementoutcomes_dump_load(backend, agnostic_load):
 
 @pytest.mark.parametrize("agnostic_load", [False, True])
 def test_circuitresult_dump_load(backend, agnostic_load):
-    c = models.Circuit(2, density_matrix=True)
+    c = Circuit(2, density_matrix=True)
     c.add(gates.M(1, 0, basis=[gates.Z, gates.X]))
     # trigger repeated execution to build the CircuitResult object
     # from samples and recover the same exact frequencies