diff --git a/src/qibo/hamiltonians/abstract.py b/src/qibo/hamiltonians/abstract.py index 86f7d559b1..966f943659 100644 --- a/src/qibo/hamiltonians/abstract.py +++ b/src/qibo/hamiltonians/abstract.py @@ -83,6 +83,21 @@ def expectation(self, state, normalize=False): # pragma: no cover """ raise_error(NotImplementedError) + @abstractmethod + def expectation_from_samples(self, freq, qubit_map=None): # pragma: no cover + """Computes the real expectation value of a diagonal observable given the frequencies when measuring in the computational basis. + + Args: + freq (collections.Counter): the keys are the observed values in binary form + and the values the corresponding frequencies, that is the number + of times each measured value/bitstring appears. + qubit_map (tuple): Mapping between frequencies and qubits. If None, [1,...,len(key)] + + Returns: + Real number corresponding to the expectation value. + """ + raise_error(NotImplementedError) + @abstractmethod def __add__(self, o): # pragma: no cover """Add operator.""" diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py index 3ef18134ae..64537c9ab4 100644 --- a/src/qibo/hamiltonians/hamiltonians.py +++ b/src/qibo/hamiltonians/hamiltonians.py @@ -2,6 +2,7 @@ from qibo.config import EINSUM_CHARS, log, raise_error from qibo.hamiltonians.abstract import AbstractHamiltonian +from qibo.symbols import Z class Hamiltonian(AbstractHamiltonian): @@ -133,6 +134,25 @@ def expectation(self, state, normalize=False): "".format(type(state)), ) + 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.") + keys = list(freq.keys()) + if qubit_map is None: + qubit_map = list(range(int(np.log2(len(obs))))) + counts = np.array(list(freq.values())) / sum(freq.values()) + expval = 0 + kl = len(qubit_map) + for j, k in enumerate(keys): + index = 0 + for i in qubit_map: + index += int(k[qubit_map.index(i)]) * 2 ** (kl - 1 - i) + expval += obs[index, index] * counts[j] + return expval + def eye(self, n=None): if n is None: n = int(self.matrix.shape[0]) @@ -534,6 +554,42 @@ def calculate_dense(self): 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: + for factor in term.factors: + if isinstance(factor, Z) == False: + raise_error( + NotImplementedError, "Observable is not a Z Pauli string." + ) + if len(term.factors) != len(set(term.factors)): + raise_error(NotImplementedError, "Z^k is not implemented since Z^2=I.") + keys = list(freq.keys()) + counts = np.array(list(freq.values())) / sum(freq.values()) + coeff = list(self.form.as_coefficients_dict().values()) + qubits = [] + for term in terms: + qubits_term = [] + for k in term.target_qubits: + qubits_term.append(k) + qubits.append(qubits_term) + if qubit_map is None: + qubit_map = list(range(len(keys[0]))) + expval = 0 + for j, q in enumerate(qubits): + subk = [] + expval_q = 0 + for i, k in enumerate(keys): + subk = [int(k[qubit_map.index(s)]) for s in q] + expval_k = 1 + if subk.count(1) % 2 == 1: + expval_k = -1 + expval_q += expval_k * counts[i] + expval += expval_q * float(coeff[j]) + return expval + def __add__(self, o): if isinstance(o, self.__class__): if self.nqubits != o.nqubits: diff --git a/src/qibo/states.py b/src/qibo/states.py index eba06884ef..972d7f6b4e 100644 --- a/src/qibo/states.py +++ b/src/qibo/states.py @@ -231,3 +231,16 @@ def apply_bitflips(self, p0, p1=None): ) noiseless_samples = self.samples() return self.backend.apply_bitflips(noiseless_samples, probs) + + def expectation_from_samples(self, observable): + """Computes the real expectation value of a diagonal observable from frequencies. + + Args: + observable (Hamiltonian/SymbolicHamiltonian): diagonal observable in the computational basis. + + Returns: + Real number corresponding to the expectation value. + """ + freq = self.frequencies(binary=True) + qubit_map = self.circuit.measurement_gate.qubits + return observable.expectation_from_samples(freq, qubit_map) diff --git a/src/qibo/tests/test_hamiltonians.py b/src/qibo/tests/test_hamiltonians.py index 759824773e..63463c2ce2 100644 --- a/src/qibo/tests/test_hamiltonians.py +++ b/src/qibo/tests/test_hamiltonians.py @@ -2,7 +2,9 @@ import numpy as np import pytest -from qibo import hamiltonians +from qibo import gates, hamiltonians +from qibo.models import Circuit +from qibo.symbols import I, Z from qibo.tests.utils import random_complex, random_sparse_matrix @@ -245,6 +247,37 @@ def test_hamiltonian_expectation_errors(backend): h.expectation("test") +def test_hamiltonian_expectation_from_samples(backend): + """Test Hamiltonian expectation value calculation.""" + obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2) + obs1 = 2 * Z(0) * Z(1) + Z(0) * Z(2) * I(3) + h0 = hamiltonians.SymbolicHamiltonian(obs0, backend=backend) + h1 = hamiltonians.SymbolicHamiltonian(obs1, backend=backend) + matrix = backend.to_numpy(h0.matrix) + c = Circuit(4) + c.add(gates.RX(0, np.random.rand())) + c.add(gates.RX(1, np.random.rand())) + c.add(gates.RX(2, np.random.rand())) + c.add(gates.RX(3, np.random.rand())) + c.add(gates.M(0, 1, 2, 3)) + nshots = 10**5 + result = c(nshots=nshots) + freq = result.frequencies(binary=True) + Obs0 = hamiltonians.Hamiltonian(3, matrix).expectation_from_samples( + freq, qubit_map=None + ) + Obs1 = h1.expectation(result.state()) + + backend.assert_allclose(Obs0, Obs1, atol=10 / np.sqrt(nshots)) + + +def test_hamiltonian_expectation_from_samples_errors(backend): + obs = random_complex((4, 4)) + h = hamiltonians.Hamiltonian(2, obs, backend=backend) + with pytest.raises(NotImplementedError): + h.expectation_from_samples(None, qubit_map=None) + + @pytest.mark.parametrize("dtype", [np.complex128, np.complex64]) @pytest.mark.parametrize( "sparse_type,dense", diff --git a/src/qibo/tests/test_hamiltonians_symbolic.py b/src/qibo/tests/test_hamiltonians_symbolic.py index 2f39f0e47a..f5e89b1b14 100644 --- a/src/qibo/tests/test_hamiltonians_symbolic.py +++ b/src/qibo/tests/test_hamiltonians_symbolic.py @@ -3,7 +3,9 @@ import pytest import sympy -from qibo import hamiltonians +from qibo import gates, hamiltonians +from qibo.models import Circuit +from qibo.symbols import I, Y, Z from qibo.tests.utils import random_complex @@ -269,6 +271,36 @@ def test_symbolic_hamiltonian_state_ev( backend.assert_allclose(local_ev, target_ev) +def test_hamiltonian_expectation_from_samples(backend): + """Test Hamiltonian expectation value calculation.""" + obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2) + obs1 = 2 * Z(0) * Z(1) + Z(0) * Z(2) * I(3) + h0 = hamiltonians.SymbolicHamiltonian(obs0, backend=backend) + h1 = hamiltonians.SymbolicHamiltonian(obs1, backend=backend) + c = Circuit(4) + c.add(gates.RX(0, np.random.rand())) + c.add(gates.RX(1, np.random.rand())) + c.add(gates.RX(2, np.random.rand())) + c.add(gates.RX(3, np.random.rand())) + c.add(gates.M(0, 1, 2, 3)) + nshots = 10**5 + result = 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)) + + +def test_hamiltonian_expectation_from_samples_errors(backend): + obs = [Z(0) * Y(1), Z(0) * Z(1) ** 3] + h1 = hamiltonians.SymbolicHamiltonian(obs[0], backend=backend) + h2 = hamiltonians.SymbolicHamiltonian(obs[1], backend=backend) + with pytest.raises(NotImplementedError): + h1.expectation_from_samples(None, qubit_map=None) + with pytest.raises(NotImplementedError): + h2.expectation_from_samples(None, qubit_map=None) + + @pytest.mark.parametrize("density_matrix", [False, True]) @pytest.mark.parametrize("calcterms", [False, True]) def test_symbolic_hamiltonian_abstract_symbol_ev(backend, density_matrix, calcterms): diff --git a/src/qibo/tests/test_states.py b/src/qibo/tests/test_states.py index f129d656f7..d46a944b73 100644 --- a/src/qibo/tests/test_states.py +++ b/src/qibo/tests/test_states.py @@ -1,7 +1,9 @@ +import numpy as np import pytest -from qibo import gates +from qibo import gates, hamiltonians from qibo.models import Circuit +from qibo.symbols import I, Z @pytest.mark.parametrize("target", range(5)) @@ -52,3 +54,26 @@ def test_state_representation_max_terms(backend, density_matrix): result.symbolic(max_terms=5) == "(0.17678+0j)|00000> + (0.17678+0j)|00001> + (0.17678+0j)|00010> + (0.17678+0j)|00011> + (0.17678+0j)|00100> + ..." ) + + +@pytest.mark.parametrize("density_matrix", [False, True]) +def test_expectation_from_samples(backend, density_matrix): + + obs0 = 2 * Z(0) * Z(1) + Z(0) * Z(2) + obs1 = 2 * Z(0) * Z(1) + Z(0) * Z(2) * I(3) + h_sym = hamiltonians.SymbolicHamiltonian(obs0, backend=backend) + h_dense = hamiltonians.Hamiltonian(3, h_sym.matrix, backend=backend) + h1 = hamiltonians.SymbolicHamiltonian(obs1, backend=backend) + c = Circuit(4) + c.add(gates.RX(0, np.random.rand())) + c.add(gates.RX(1, np.random.rand())) + c.add(gates.RX(2, np.random.rand())) + c.add(gates.RX(3, np.random.rand())) + c.add(gates.M(0, 1, 2)) + nshots = 10**5 + result = c(nshots=nshots) + expval_sym = result.expectation_from_samples(h_sym) + expval_dense = result.expectation_from_samples(h_dense) + expval = h1.expectation(result.state()) + backend.assert_allclose(expval_sym, expval_dense) + backend.assert_allclose(expval_sym, expval, atol=10 / np.sqrt(nshots))