From 8abe31fecd5aa44e1611cd69e5bad325c56106ce Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Sat, 24 Jun 2023 11:29:54 +0300
Subject: [PATCH 1/7] Add nqubits and errors in SymbolicHamiltonian

---
 src/qibo/hamiltonians/hamiltonians.py | 34 +++++++--------
 tests/test_hamiltonians_symbolic.py   | 59 +++++++++++++++++++++++----
 2 files changed, 69 insertions(+), 24 deletions(-)

diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py
index 804f9f261d..41522d310d 100644
--- a/src/qibo/hamiltonians/hamiltonians.py
+++ b/src/qibo/hamiltonians/hamiltonians.py
@@ -1,3 +1,6 @@
+from itertools import chain
+
+import numpy as np
 import sympy
 
 from qibo.config import EINSUM_CHARS, log, raise_error
@@ -135,8 +138,6 @@ def expectation(self, state, normalize=False):
             )
 
     def expectation_from_samples(self, freq, qubit_map=None):
-        import numpy as np
-
         obs = self.matrix
         if np.count_nonzero(obs - np.diag(np.diagonal(obs))) != 0:
             raise_error(NotImplementedError, "Observable is not diagonal.")
@@ -318,7 +319,7 @@ class SymbolicHamiltonian(AbstractHamiltonian):
             :meth:`qibo.hamiltonians.models.MaxCut` Hamiltonian.
     """
 
-    def __init__(self, form=None, symbol_map={}, backend=None):
+    def __init__(self, form=None, nqubits=None, symbol_map={}, backend=None):
         super().__init__()
         self._form = None
         self._terms = None
@@ -339,6 +340,8 @@ def __init__(self, form=None, symbol_map={}, backend=None):
             self.backend = backend
         if form is not None:
             self.form = form
+        if nqubits is not None:
+            self.nqubits = nqubits
 
     @property
     def dense(self):
@@ -409,9 +412,9 @@ def terms(self):
                     terms.append(term)
                 else:
                     self.constant += term.coefficient
-            assert (
-                self.nqubits == max(q for term in terms for q in term.target_qubits) + 1
-            )
+            # assert (
+            #    self.nqubits == max(q for term in terms for q in term.target_qubits) + 1
+            # )
             self._terms = terms
         return self._terms
 
@@ -454,8 +457,6 @@ def _get_symbol_matrix(self, term):
             Numerical matrix corresponding to the given expression as a numpy
             array of size ``(2 ** self.nqubits, 2 ** self.nqubits).
         """
-        from numpy import eye
-
         if isinstance(term, sympy.Add):
             # symbolic op for addition
             result = sum(
@@ -492,7 +493,7 @@ def _get_symbol_matrix(self, term):
                 if not isinstance(matrix, self.backend.tensor_types):
                     # symbols that do not correspond to quantum operators
                     # for example parameters in the MaxCut Hamiltonian
-                    result = complex(matrix) * eye(2**self.nqubits)
+                    result = complex(matrix) * np.eye(2**self.nqubits)
                 else:
                     # if we do not have a Qibo symbol we construct one and use
                     # :meth:`qibo.core.terms.SymbolicTerm.full_matrix`.
@@ -502,7 +503,7 @@ def _get_symbol_matrix(self, term):
             # if the term is number we should return in the form of identity
             # matrix because in expressions like `1 + Z`, `1` is not correspond
             # to the float 1 but the identity operator (matrix)
-            result = complex(term) * eye(2**self.nqubits)
+            result = complex(term) * np.eye(2**self.nqubits)
 
         else:
             raise_error(
@@ -522,10 +523,6 @@ def _calculate_dense_from_form(self):
 
     def _calculate_dense_from_terms(self):
         """Calculates equivalent :class:`qibo.core.hamiltonians.Hamiltonian` using the term representation."""
-        from itertools import chain
-
-        import numpy as np
-
         if 2 * self.nqubits > len(EINSUM_CHARS):  # pragma: no cover
             # case not tested because it only happens in large examples
             raise_error(NotImplementedError, "Not enough einsum characters.")
@@ -555,8 +552,6 @@ def expectation(self, state, normalize=False):
         return Hamiltonian.expectation(self, state, normalize)
 
     def expectation_from_samples(self, freq, qubit_map=None):
-        import numpy as np
-
         terms = self.terms
         for term in terms:
             # pylint: disable=E1101
@@ -745,6 +740,13 @@ def __matmul__(self, o):
 
         if isinstance(o, self.backend.tensor_types):
             rank = len(tuple(o.shape))
+            state_qubits = int(np.log2(int(o.shape[0])))
+            if state_qubits != self.nqubits:
+                raise_error(
+                    ValueError,
+                    f"Cannot multiply Hamiltonian on {self.nqubits} qubits to "
+                    f"state of {state_qubits} qubits.",
+                )
             if rank == 1:  # state vector
                 return self.apply_gates(o)
             elif rank == 2:  # density matrix
diff --git a/tests/test_hamiltonians_symbolic.py b/tests/test_hamiltonians_symbolic.py
index c9d208c047..757c5a9705 100644
--- a/tests/test_hamiltonians_symbolic.py
+++ b/tests/test_hamiltonians_symbolic.py
@@ -35,9 +35,13 @@ def test_symbolic_hamiltonian_errors(backend):
     Z, X = sympy.Symbol("Z"), sympy.Symbol("X")
     symbol_map = {Z: (0, matrices.Z)}
     with pytest.raises(ValueError):
-        ham = hamiltonians.SymbolicHamiltonian(Z * X, symbol_map, backend=backend)
+        ham = hamiltonians.SymbolicHamiltonian(
+            Z * X, symbol_map=symbol_map, backend=backend
+        )
     # Invalid operation in Hamiltonian expresion
-    ham = hamiltonians.SymbolicHamiltonian(sympy.cos(Z), symbol_map, backend=backend)
+    ham = hamiltonians.SymbolicHamiltonian(
+        sympy.cos(Z), symbol_map=symbol_map, backend=backend
+    )
     with pytest.raises(TypeError):
         dense = ham.dense
 
@@ -248,7 +252,7 @@ def test_symbolic_hamiltonian_matmul(backend, nqubits, density_matrix, calcterms
 @pytest.mark.parametrize("nqubits,normalize", [(3, False), (4, False)])
 @pytest.mark.parametrize("calcterms", [False, True])
 @pytest.mark.parametrize("calcdense", [False, True])
-def test_symbolic_hamiltonian_state_ev(
+def test_symbolic_hamiltonian_state_expectation(
     backend, nqubits, normalize, calcterms, calcdense
 ):
     local_ham = (
@@ -272,6 +276,45 @@ def test_symbolic_hamiltonian_state_ev(
     backend.assert_allclose(local_ev, target_ev)
 
 
+@pytest.mark.parametrize("give_nqubits", [False, True])
+@pytest.mark.parametrize("calcterms", [False, True])
+@pytest.mark.parametrize("calcdense", [False, True])
+def test_symbolic_hamiltonian_state_expectation_different_nqubits(
+    backend, give_nqubits, calcterms, calcdense
+):
+    expr = symbolic_tfim(3, h=1.0)
+    if give_nqubits:
+        local_ham = hamiltonians.SymbolicHamiltonian(expr, nqubits=5, backend=backend)
+    else:
+        local_ham = hamiltonians.SymbolicHamiltonian(expr, backend=backend)
+    if calcterms:
+        _ = local_ham.terms
+    if calcdense:
+        _ = local_ham.dense
+
+    dense_ham = hamiltonians.TFIM(3, h=1.0, backend=backend)
+    dense_matrix = np.kron(backend.to_numpy(dense_ham.matrix), np.eye(4))
+    dense_ham = hamiltonians.Hamiltonian(5, dense_matrix, backend=backend)
+
+    if give_nqubits:
+        state = backend.cast(random_complex((2**5,)))
+        local_ev = local_ham.expectation(state)
+        target_ev = dense_ham.expectation(state)
+        backend.assert_allclose(local_ev, target_ev)
+
+        state = random_complex((2**5,))
+        local_ev = local_ham.expectation(state)
+        target_ev = dense_ham.expectation(state)
+        backend.assert_allclose(local_ev, target_ev)
+    else:
+        state = backend.cast(random_complex((2**5,)))
+        with pytest.raises(ValueError):
+            local_ev = local_ham.expectation(state)
+        state = random_complex((2**5,))
+        with pytest.raises(ValueError):
+            local_ev = local_ham.expectation(state)
+
+
 def test_hamiltonian_expectation_from_samples(backend):
     """Test Hamiltonian expectation value calculation."""
     obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2)
@@ -285,11 +328,11 @@ def test_hamiltonian_expectation_from_samples(backend):
     c.add(gates.RX(3, np.random.rand()))
     c.add(gates.M(0, 1, 2, 3))
     nshots = 10**5
-    result = c(nshots=nshots)
+    result = backend.execute_circuit(c, nshots=nshots)
     freq = result.frequencies(binary=True)
-    Obs0 = h0.expectation_from_samples(freq, qubit_map=None)
-    Obs1 = h1.expectation(result.state())
-    backend.assert_allclose(Obs0, Obs1, atol=10 / np.sqrt(nshots))
+    ev0 = h0.expectation_from_samples(freq, qubit_map=None)
+    ev1 = h1.expectation(result.state())
+    backend.assert_allclose(ev0, ev1, atol=20 / np.sqrt(nshots))
 
 
 def test_hamiltonian_expectation_from_samples_errors(backend):
@@ -341,7 +384,7 @@ def test_trotter_hamiltonian_operation_errors(backend):
         h = h1 * "test"
     with pytest.raises(NotImplementedError):
         h = h1 @ "test"
-    with pytest.raises(NotImplementedError):
+    with pytest.raises(ValueError):
         h = h1 @ np.ones((2, 2, 2, 2))
     h2 = hamiltonians.XXZ(3, dense=False, backend=backend)
     with pytest.raises(NotImplementedError):

From ed75760c15d134fa5c998add199b673ec51b98dc Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Sat, 24 Jun 2023 11:51:34 +0300
Subject: [PATCH 2/7] Fix symbol map

---
 src/qibo/hamiltonians/hamiltonians.py | 4 +++-
 src/qibo/hamiltonians/models.py       | 2 +-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py
index 41522d310d..e090dfd7e9 100644
--- a/src/qibo/hamiltonians/hamiltonians.py
+++ b/src/qibo/hamiltonians/hamiltonians.py
@@ -91,7 +91,9 @@ def from_symbolic(cls, symbolic_hamiltonian, symbol_map, backend=None):
             "deprecated. Please use `SymbolicHamiltonian` and Qibo symbols "
             "to construct Hamiltonians using symbols."
         )
-        return SymbolicHamiltonian(symbolic_hamiltonian, symbol_map, backend=backend)
+        return SymbolicHamiltonian(
+            symbolic_hamiltonian, symbol_map=symbol_map, backend=backend
+        )
 
     def eigenvalues(self, k=6):
         if self._eigenvalues is None:
diff --git a/src/qibo/hamiltonians/models.py b/src/qibo/hamiltonians/models.py
index f1171e27f5..1fa030fac6 100644
--- a/src/qibo/hamiltonians/models.py
+++ b/src/qibo/hamiltonians/models.py
@@ -191,7 +191,7 @@ def MaxCut(nqubits, dense=True, backend=None):
     smap = {s: (i, matrices.Z) for i, s in enumerate(Z)}
     smap.update({s: (i, v[i]) for i, s in enumerate(V)})
 
-    ham = SymbolicHamiltonian(sham, smap, backend=backend)
+    ham = SymbolicHamiltonian(sham, symbol_map=smap, backend=backend)
     if dense:
         return ham.dense
     return ham

From f23e2c24426c89fb9b78100472a8965985bb8e98 Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Sat, 24 Jun 2023 12:28:14 +0300
Subject: [PATCH 3/7] Add example for expectation from samples

---
 doc/source/api-reference/qibo.rst             |  2 -
 doc/source/code-examples/advancedexamples.rst | 71 +++++++++++++++++++
 2 files changed, 71 insertions(+), 2 deletions(-)

diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index c89102799f..8d1178aa6d 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -851,7 +851,6 @@ values using a sparse Hamiltonian matrix.
 .. autoclass:: qibo.hamiltonians.Hamiltonian
     :members:
     :member-order: bysource
-    :noindex:
 
 
 Symbolic Hamiltonian
@@ -866,7 +865,6 @@ For more information on constructing Hamiltonians using symbols we refer to the
 .. autoclass:: qibo.hamiltonians.SymbolicHamiltonian
     :members:
     :member-order: bysource
-    :noindex:
 
 
 When a :class:`qibo.hamiltonians.SymbolicHamiltonian` is used for time
diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 82c95771a1..fd6e7674cb 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1897,3 +1897,74 @@ constructing each symbol:
 
     form = Z(0, commutative=True) * Z(1, commutative=True) + Z(1, commutative=True) * Z(2, commutative=True)
     ham = hamiltonians.SymbolicHamiltonian(form)
+
+
+.. _hamexpectation-example:
+
+How to calculate expectation values using samples?
+--------------------------------------------------
+
+It is possible to calculate the expectation value of a :class:`qibo.hamiltonians.Hamiltonian`
+on a given state using the :meth:`qibo.hamiltonians.Hamiltonian.expectation` method,
+which can be called on a state or density matrix. For example
+
+
+.. testcode::
+
+    import numpy as np
+    from qibo.models import Circuit
+    from qibo.hamiltonians import XXZ
+
+    circuit = Circuit(4)
+    circuit.add(gates.H(i) for i in range(4))
+    c.add(gates.CNOT(0, 1))
+    c.add(gates.CNOT(1, 2))
+    c.add(gates.CNOT(2, 3))
+
+    hamiltonian = XXZ(4)
+
+    result = circuit()
+    expectation_value = hamiltonian.expectation(result.state())
+
+In this example, the circuit will be simulated to obtain the final state vector
+and the corresponding expectation value will be calculated through exact matrix
+multiplication with the Hamiltonian matrix.
+If a :class:`qibo.hamiltonians.SymbolicHamiltonian` is used instead, the expectation
+value will be calculated as a sum of expectation values of local terms, allowing
+calculations of more qubits with lower memory consumption. The calculation of each
+local term still requires the state vector.
+
+When executing a circuit on real hardware, usually only measurements of the state are
+available, not the state vector. Qibo provides :meth:`qibo.hamiltonians.Hamiltonian.expectation_from_samples`
+to allow calculation of expectation values directly from such samples:
+
+
+.. testcode::
+
+    import numpy as np
+    from qibo.models import Circuit
+    from qibo.hamiltonians import XXZ
+
+    circuit = Circuit(4)
+    circuit.add(gates.H(i) for i in range(4))
+    c.add(gates.CNOT(0, 1))
+    c.add(gates.CNOT(1, 2))
+    c.add(gates.CNOT(2, 3))
+
+    hamiltonian = XXZ(4)
+
+    result = circuit()
+    expectation_value = hamiltonian.expectation_from_samples(result.frequencies())
+
+
+This example simulates the circuit similarly to the previous one but calculates
+the expectation value using the frequencies of shots, instead of the exact state vector.
+This can also be invoked directly from the ``result`` object:
+
+.. testcode::
+
+    expectation_value = result.expectation_from_samples(hamiltonian)
+
+
+The expectation from samples currently works only for Hamiltonians that are diagonal in
+the computational basis.

From f3e3b7dacd413d4da99fab2c9e9ae0d891f50749 Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Sat, 24 Jun 2023 12:53:22 +0300
Subject: [PATCH 4/7] Fix coverage

---
 src/qibo/hamiltonians/hamiltonians.py | 15 ++++++---------
 tests/test_hamiltonians_symbolic.py   |  2 +-
 2 files changed, 7 insertions(+), 10 deletions(-)

diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py
index e090dfd7e9..351a39f559 100644
--- a/src/qibo/hamiltonians/hamiltonians.py
+++ b/src/qibo/hamiltonians/hamiltonians.py
@@ -414,9 +414,6 @@ def terms(self):
                     terms.append(term)
                 else:
                     self.constant += term.coefficient
-            # assert (
-            #    self.nqubits == max(q for term in terms for q in term.target_qubits) + 1
-            # )
             self._terms = terms
         return self._terms
 
@@ -742,6 +739,11 @@ def __matmul__(self, o):
 
         if isinstance(o, self.backend.tensor_types):
             rank = len(tuple(o.shape))
+            if rank not in (1, 2):
+                raise_error(
+                    NotImplementedError,
+                    "Cannot multiply Hamiltonian with " "rank-{} tensor.".format(rank),
+                )
             state_qubits = int(np.log2(int(o.shape[0])))
             if state_qubits != self.nqubits:
                 raise_error(
@@ -751,13 +753,8 @@ def __matmul__(self, o):
                 )
             if rank == 1:  # state vector
                 return self.apply_gates(o)
-            elif rank == 2:  # density matrix
+            else:  # density matrix
                 return self.apply_gates(o, density_matrix=True)
-            else:
-                raise_error(
-                    NotImplementedError,
-                    "Cannot multiply Hamiltonian with " "rank-{} tensor.".format(rank),
-                )
 
         raise_error(
             NotImplementedError,
diff --git a/tests/test_hamiltonians_symbolic.py b/tests/test_hamiltonians_symbolic.py
index 757c5a9705..dea02fc124 100644
--- a/tests/test_hamiltonians_symbolic.py
+++ b/tests/test_hamiltonians_symbolic.py
@@ -384,7 +384,7 @@ def test_trotter_hamiltonian_operation_errors(backend):
         h = h1 * "test"
     with pytest.raises(NotImplementedError):
         h = h1 @ "test"
-    with pytest.raises(ValueError):
+    with pytest.raises(NotImplementedError):
         h = h1 @ np.ones((2, 2, 2, 2))
     h2 = hamiltonians.XXZ(3, dense=False, backend=backend)
     with pytest.raises(NotImplementedError):

From 0355dc0eff42f134d8608ab0f59715b86e708c5f Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Sat, 24 Jun 2023 13:24:47 +0300
Subject: [PATCH 5/7] Fix nshots

Co-authored-by: Sergi Ramos <36544268+igres26@users.noreply.github.com>
---
 doc/source/code-examples/advancedexamples.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index fd6e7674cb..192b3e2d45 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1953,7 +1953,7 @@ to allow calculation of expectation values directly from such samples:
 
     hamiltonian = XXZ(4)
 
-    result = circuit()
+    result = circuit(nshots=1024)
     expectation_value = hamiltonian.expectation_from_samples(result.frequencies())
 
 

From af2560674ea8016144ff2323faac58fd56216642 Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Mon, 26 Jun 2023 12:38:48 +0300
Subject: [PATCH 6/7] Fix doc code examples

---
 doc/source/code-examples/advancedexamples.rst | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index 192b3e2d45..c066e32ebd 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1911,7 +1911,7 @@ which can be called on a state or density matrix. For example
 
 .. testcode::
 
-    import numpy as np
+    from qibo import gates
     from qibo.models import Circuit
     from qibo.hamiltonians import XXZ
 
@@ -1941,9 +1941,9 @@ to allow calculation of expectation values directly from such samples:
 
 .. testcode::
 
-    import numpy as np
+    from qibo import gates
     from qibo.models import Circuit
-    from qibo.hamiltonians import XXZ
+    from qibo.hamiltonians import Z
 
     circuit = Circuit(4)
     circuit.add(gates.H(i) for i in range(4))
@@ -1951,7 +1951,7 @@ to allow calculation of expectation values directly from such samples:
     c.add(gates.CNOT(1, 2))
     c.add(gates.CNOT(2, 3))
 
-    hamiltonian = XXZ(4)
+    hamiltonian = Z(4)
 
     result = circuit(nshots=1024)
     expectation_value = hamiltonian.expectation_from_samples(result.frequencies())

From 23ce1011701d05ca93f94fe5980ea5b5103c0622 Mon Sep 17 00:00:00 2001
From: Stavros Efthymiou <35475381+stavros11@users.noreply.github.com>
Date: Wed, 28 Jun 2023 08:51:00 +0300
Subject: [PATCH 7/7] Doc example fixes

---
 doc/source/code-examples/advancedexamples.rst | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst
index c066e32ebd..81c712f910 100644
--- a/doc/source/code-examples/advancedexamples.rst
+++ b/doc/source/code-examples/advancedexamples.rst
@@ -1917,9 +1917,9 @@ which can be called on a state or density matrix. For example
 
     circuit = Circuit(4)
     circuit.add(gates.H(i) for i in range(4))
-    c.add(gates.CNOT(0, 1))
-    c.add(gates.CNOT(1, 2))
-    c.add(gates.CNOT(2, 3))
+    circuit.add(gates.CNOT(0, 1))
+    circuit.add(gates.CNOT(1, 2))
+    circuit.add(gates.CNOT(2, 3))
 
     hamiltonian = XXZ(4)
 
@@ -1947,9 +1947,10 @@ to allow calculation of expectation values directly from such samples:
 
     circuit = Circuit(4)
     circuit.add(gates.H(i) for i in range(4))
-    c.add(gates.CNOT(0, 1))
-    c.add(gates.CNOT(1, 2))
-    c.add(gates.CNOT(2, 3))
+    circuit.add(gates.CNOT(0, 1))
+    circuit.add(gates.CNOT(1, 2))
+    circuit.add(gates.CNOT(2, 3))
+    circuit.add(gates.M(*range(4)))
 
     hamiltonian = Z(4)