Adiabatic evolution class state callback returns double the expected states
When adding state callbacks to the adiabatic evolution function the resulting state object has repeated states for each state.

To reproduce one can follow the code depicted in the images which also shows the issue.

![Screenshot 2023-07-19 at 10 35 22](https://github.com/qiboteam/qibo/assets/139958400/75d2bdec-a402-4414-b7d5-c77cad580a37)
![Screenshot 2023-07-19 at 10 37 07](https://github.com/qiboteam/qibo/assets/139958400/6a81ff49-3394-4465-9f57-0fd41f0f71e0)

Im using the latest qibo package (v 0.1.14) . I would expect both lengths to match or in this case to differ by 1 if counting the initial state when t=0 which it does.


diff --git a/src/qibo/callbacks.py b/src/qibo/callbacks.py
index aa310fe3..f68e8d23 100644
--- a/src/qibo/callbacks.py
+++ b/src/qibo/callbacks.py
@@ -105,8 +105,6 @@ class EntanglementEntropy(Callback):
 
     @Callback.nqubits.setter
     def nqubits(self, n: int):
-        from qibo import gates
-
         if self._nqubits is not None and self._nqubits != n:
             raise_error(
                 RuntimeError,
@@ -115,9 +113,9 @@ class EntanglementEntropy(Callback):
         self._nqubits = n
         if self.partition is None:
             self.partition = list(range(n // 2 + n % 2))
-        if len(self.partition) <= self.nqubits // 2:
+        if len(self.partition) <= self._nqubits // 2:
             self.partition = [
-                i for i in range(self.nqubits) if i not in set(self.partition)
+                i for i in range(self._nqubits) if i not in set(self.partition)
             ]
 
     def apply(self, backend, state):
@@ -175,10 +173,14 @@ class Norm(Callback):
     """
 
     def apply(self, backend, state):
-        return backend.calculate_norm(state)
+        norm = backend.calculate_norm(state)
+        self.append(norm)
+        return norm
 
     def apply_density_matrix(self, backend, state):
-        return backend.calculate_norm_density_matrix(state)
+        norm = backend.calculate_norm_density_matrix(state)
+        self.append(norm)
+        return norm
 
 
 class Overlap(Callback):
@@ -200,10 +202,14 @@ class Overlap(Callback):
         self.state = state
 
     def apply(self, backend, state):
-        return backend.calculate_overlap(self.state, state)
+        overlap = backend.calculate_overlap(self.state, state)
+        self.append(overlap)
+        return overlap
 
-    def apply_density_matrix(self, backend, state):
-        return backend.calculate_overlap_density_matrix(self.state, state)
+    def apply_density_matrix(self, backend, state):  # pragma: no cover
+        overlap = backend.calculate_overlap_density_matrix(self.state, state)
+        self.append(overlap)
+        return overlap
 
 
 class Energy(Callback):
@@ -229,11 +235,15 @@ class Energy(Callback):
 
     def apply(self, backend, state):
         assert type(self.hamiltonian.backend) == type(backend)
-        return self.hamiltonian.expectation(state)
+        expectation = self.hamiltonian.expectation(state)
+        self.append(expectation)
+        return expectation
 
     def apply_density_matrix(self, backend, state):
         assert type(self.hamiltonian.backend) == type(backend)
-        return self.hamiltonian.expectation(state)
+        expectation = self.hamiltonian.expectation(state)
+        self.append(expectation)
+        return expectation
 
 
 class Gap(Callback):
@@ -314,12 +324,16 @@ class Gap(Callback):
         hamiltonian.eigenvectors()
         eigvals = hamiltonian.eigenvalues()
         if isinstance(self.mode, int):
-            return backend.np.real(eigvals[self.mode])
+            gap = backend.np.real(eigvals[self.mode])
+            self.append(gap)
+            return gap
 
         # case: self.mode == "gap"
         excited = 1
         gap = backend.np.real(eigvals[excited] - eigvals[0])
+
         if not self.check_degenerate:
+            self.append(gap)
             return gap
 
         while backend.np.less(gap, EIGVAL_CUTOFF):
@@ -330,6 +344,7 @@ class Gap(Callback):
                 "The Hamiltonian is degenerate. Using eigenvalue {} "
                 "to calculate gap.".format(excited)
             )
+        self.append(gap)
         return gap
 
     def apply_density_matrix(self, backend, state):
diff --git a/src/qibo/models/evolution.py b/src/qibo/models/evolution.py
index 42c43be2..9a9cee88 100644
--- a/src/qibo/models/evolution.py
+++ b/src/qibo/models/evolution.py
@@ -89,7 +89,10 @@ class StateEvolution:
     def _create_calculate_callbacks(self, accelerators):
         def calculate_callbacks(state):
             for callback in self.callbacks:
-                callback.append(callback.apply(self.backend, state))
+                callback.nqubits = self.nqubits
+                # by executing callbacks.apply we also append the object to history
+                # see callbacks module for this
+                callback.apply(self.backend, state)
 
         if accelerators is None:
             return calculate_callbacks
['src/qibo/models/evolution.py', 'src/qibo/callbacks.py']