feat: separated the two expectation functions

doc/source/code-examples/advancedexamples.rst

@@ -2058,7 +2058,7 @@ diagonalised simultaneously, one has to calculate the expectation value starting
    # while calculating the expectation value
 
    hamiltonian = SymbolicHamiltonian(3 * Z(2) * (1 - X(1)) ** 2 - (Y(0) * X(3)) / 2, nqubits=4)
-   expectation_value = hamiltonian.expectation_from_samples(circuit)
+   expectation_value = hamiltonian.expectation_from_circuit(circuit)
 
 What is happening under the hood in this case, is that the expectation value is calculated for each term
 individually by measuring the circuit in the correct (rotated) basis. All the contributions are then

src/qibo/hamiltonians/hamiltonians.py

@@ -2,7 +2,7 @@
 
 from itertools import chain
 from math import prod
-from typing import Optional, Union
+from typing import Optional
 
 import numpy as np
 import sympy
@@ -556,26 +556,34 @@ def calculate_dense(self):
     def expectation(self, state, normalize=False):
         return Hamiltonian.expectation(self, state, normalize)
 
-    def _exp_from_circuit(
-        self, circuit: "Circuit", qubit_map: dict, nshots: int = None
-    ):
+    def expectation_from_circuit(
+        self, circuit: "Circuit", qubit_map: dict = None, nshots: int = 1000
+    ) -> float:
         """
         Calculate the expectation value from a circuit.
-        This allows for non diagonal observables. Each term of the observable is
-        treated separately, by measuring in the correct basis and re-executing the
-        circuit.
+        This even works for observables not completely diagonal in the computational
+        basis, but only diagonal at a term level in a defined basis. Namely, for
+        an observable of the form :math:``H = \\sum_i H_i``, where each :math:``H_i``
+        consists in a `n`-qubits pauli operator :math:`P_0 \\otimes P_1 \\otimes \\cdots \\otimes P_n`,
+        the expectation value is computed by rotating the input circuit in the suitable
+        basis for each term :math:``H_i`` thus extracting the `term-wise` expectations
+        that are then summed to build the global expectation value.
+        Each term of the observable is treated separately, by measuring in the correct
+        basis and re-executing the circuit.
 
         Args:
-            circuit (Circuit): input frequencies.
-            qubit_map (dict): qubit map.
+            circuit (Circuit): input circuit.
+            qubit_map (dict): qubit map, defaults to ``{0: 0, 1: 1, ..., n: n}``
+            nshots (int): number of shots, defaults to 1000.
 
         Returns:
             (float): the calculated expectation value.
         """
         from qibo import Circuit, gates
 
-        if nshots is None:  # pragma: no cover
-            nshots = 1000
+        if qubit_map is None:
+            qubit_map = list(range(self.nqubits))
+
         rotated_circuits = []
         coefficients = []
         Z_observables = []
@@ -614,20 +622,20 @@ def _exp_from_circuit(
         ]
         return sum(
             [
-                coeff * obs._exp_from_freq(freq, qmap)
+                coeff * obs.expectation_from_samples(freq, qmap)
                 for coeff, freq, obs, qmap in zip(
                     coefficients, frequencies, Z_observables, q_maps
                 )
             ]
         )
 
-    def _exp_from_freq(self, freq: dict, qubit_map: dict) -> float:
+    def expectation_from_samples(self, freq: dict, qubit_map: dict = None) -> float:
         """
-        Calculate the expectation value from some frequencies.
+        Calculate the expectation value from the samples.
         The observable has to be diagonal in the computational basis.
 
         Args:
-            freq (dict): input frequencies.
+            freq (dict): input frequencies of the samples.
             qubit_map (dict): qubit map.
 
         Returns:
@@ -641,6 +649,9 @@ def _exp_from_freq(self, freq: dict, qubit_map: dict) -> float:
                         NotImplementedError, "Observable is not a Z Pauli string."
                     )
 
+        if qubit_map is None:
+            qubit_map = list(range(self.nqubits))
+
         keys = list(freq.keys())
         counts = self.backend.cast(list(freq.values()), self.backend.precision) / sum(
             freq.values()
@@ -660,38 +671,6 @@ def _exp_from_freq(self, freq: dict, qubit_map: dict) -> float:
         )
         return self.backend.np.sum(expvals @ counts.T) + self.constant.real
 
-    def expectation_from_samples(
-        self, data: Union[dict, "Circuit"], qubit_map: dict = None, nshots: int = None
-    ):
-        """
-        Calculate the expectation value starting from the samples. This even works
-        for observables not completely diagonal in the computational basis, but only
-        diagonal at a term level in a defined basis. Namely, for an observable of the
-        form :math:``H = \\sum_i H_i``, where each :math:``H_i`` consists in a `n`-qubits
-        pauli operator :math:`P_0 \\otimes P_1 \\otimes \\cdots \\otimes P_n`, the
-        expectation value is computed by rotating the input circuit in the suitable
-        basis for each term :math:``H_i`` thus extracting the `term-wise` expectations
-        that are then summed to build the global expectation value.
-
-        Args:
-            data (dict | :class:`qibo.models.Circuit`): either the ``dict`` of the
-            frequencies of the samples or a :class:`qibo.models.Circuit` where to
-            extract the samples from. A :class:`qibo.models.Circuit` is needed in
-            case the observable is not diagonal in the computational basis.
-            qubit_map (dict): qubit map.
-            nshots (int, optional): number of shots for the expecation value
-            calculation. This is used only when a :class:`qibo.models.Circuit`
-            is passed as input.
-
-        Returns:
-            (float) the computed expectation value.
-        """
-        if qubit_map is None:
-            qubit_map = list(range(self.nqubits))
-        if isinstance(data, dict):
-            return self._exp_from_freq(data, qubit_map)
-        return self._exp_from_circuit(data, qubit_map, nshots)
-
     def __add__(self, o):
         if isinstance(o, self.__class__):
             if self.nqubits != o.nqubits:

tests/test_hamiltonians.py

@@ -262,38 +262,40 @@ def test_hamiltonian_expectation_errors(backend):
         h.expectation("test")
 
 
-@pytest.mark.parametrize(
-    "observable",
-    [
-        2 * Z(0) * (1 - Z(1)) ** 2 + Z(0) * Z(2),
-        X(0) * Z(1) + Y(0) * X(2) / 2 - Z(0) * (1 - Y(1)) ** 3,
-    ],
-)
-def test_hamiltonian_expectation_from_samples(backend, observable):
-    """Test Hamiltonian expectation value calculation."""
-    backend.set_seed(12)
+def non_exact_expectation_test_setup(backend, observable):
 
     nqubits = 3
-    nshots = 4 * 10**6
     c = Circuit(nqubits)
     for q in range(nqubits):
         c.add(gates.RX(q, np.random.rand()))
 
     H = hamiltonians.SymbolicHamiltonian(observable, nqubits=nqubits, backend=backend)
-    matrix = H.matrix
-    diagonal = (
-        backend.np.count_nonzero(matrix - backend.np.diag(backend.np.diagonal(matrix)))
-        == 0
-    )
     final_state = backend.execute_circuit(c.copy(True)).state()
     exp = H.expectation(final_state)
-    if not diagonal:
-        exp_from_samples = H.expectation_from_samples(c, nshots=nshots)
-    else:
-        c.add(gates.M(*range(nqubits)))
-        freq = backend.execute_circuit(c, nshots=nshots).frequencies()
-        exp_from_samples = H.expectation_from_samples(freq)
+    return exp, H, c
+
+
+def test_hamiltonian_expectation_from_samples(backend):
+    """Test Hamiltonian expectation value calculation."""
+    backend.set_seed(12)
+
+    nshots = 4 * 10**6
+    observable = 2 * Z(0) * (1 - Z(1)) ** 2 + Z(0) * Z(2)
+    exp, H, c = non_exact_expectation_test_setup(backend, observable)
+    c.add(gates.M(*range(c.nqubits)))
+    freq = backend.execute_circuit(c, nshots=nshots).frequencies()
+    exp_from_samples = H.expectation_from_samples(freq)
+    backend.assert_allclose(exp, exp_from_samples, atol=1e-2)
 
+
+def test_hamiltonian_expectation_from_circuit(backend):
+    """Test Hamiltonian expectation value calculation."""
+    backend.set_seed(12)
+
+    nshots = 4 * 10**6
+    observable = X(0) * Z(1) + Y(0) * X(2) / 2 - Z(0) * (1 - Y(1)) ** 3
+    exp, H, c = non_exact_expectation_test_setup(backend, observable)
+    exp_from_samples = H.expectation_from_circuit(c, nshots=nshots)
     backend.assert_allclose(exp, exp_from_samples, atol=1e-2)