Merge pull request #59 from qiboteam/expectation-from-samples

Improvements to expectation_from_samples code

doc/source/api-reference/measurement.rst

@@ -2,6 +2,10 @@
 Measurement
 ===========
 
-This section covers the API reference for obtaining the expectation value of some (molecular) Hamiltonian
+This section covers the API reference for obtaining the expectation value of a (molecular) Hamiltonian
 
-.. autofunction:: qibochem.measurement.expectation.expectation
+.. autofunction:: qibochem.measurement.optimization.measurement_basis_rotations
+
+.. autofunction:: qibochem.measurement.optimization.allocate_shots
+
+.. autofunction:: qibochem.measurement.expectation

src/qibochem/measurement/expectation.py

@@ -1,57 +1,107 @@
+from functools import reduce
+
+import numpy as np
 import qibo
-from qibo import gates
 from qibo.hamiltonians import SymbolicHamiltonian
 from qibo.symbols import Z
 
+from qibochem.measurement.optimization import (
+    allocate_shots,
+    measurement_basis_rotations,
+)
+
+
+def symbolic_term_to_symbol(symbolic_term):
+    """Convert a single Pauli word in the form of a Qibo SymbolicTerm to a Qibo Symbol"""
+    return symbolic_term.coefficient * reduce(lambda x, y: x * y, symbolic_term.factors, 1.0)
+
+
+def pauli_term_measurement_expectation(pauli_term, frequencies):
+    """
+    Calculate the expectation value of a single general Pauli string for some measurement frequencies
+
+    Args:
+        pauli_term (SymbolicTerm): Single general Pauli term, e.g. X0*Z2
+        frequencies: Measurement frequencies, taken from MeasurementOutcomes.frequencies(binary=True)
+
+    Returns:
+        float: Expectation value of pauli_term
+    """
+    # Replace every (non-I) Symbol with Z, then include the term coefficient
+    pauli_z = [Z(int(factor.target_qubit)) for factor in pauli_term.factors if factor.name[0] != "I"]
+    z_only_ham = SymbolicHamiltonian(pauli_term.coefficient * reduce(lambda x, y: x * y, pauli_z, 1.0))
+    # Can now apply expectation_from_samples directly
+    return z_only_ham.expectation_from_samples(frequencies)
+
 
 def expectation(
-    circuit: qibo.models.Circuit, hamiltonian: SymbolicHamiltonian, from_samples=False, n_shots=1000
+    circuit: qibo.models.Circuit,
+    hamiltonian: SymbolicHamiltonian,
+    from_samples: bool = False,
+    n_shots: int = 1000,
+    group_pauli_terms=None,
+    n_shots_per_pauli_term: bool = True,
+    shot_allocation=None,
 ) -> float:
     """
     Calculate expectation value of some Hamiltonian using either the state vector or sample measurements from running a
     quantum circuit
 
     Args:
         circuit (qibo.models.Circuit): Quantum circuit ansatz
-        hamiltonian (SymbolicHamiltonian): Molecular Hamiltonian
-        from_samples (Boolean): Whether the expectation value calculation uses samples or the simulated
+        hamiltonian (qibo.hamiltonians.SymbolicHamiltonian): Molecular Hamiltonian
+        from_samples (bool): Whether the expectation value calculation uses samples or the simulated
             state vector. Default: ``False``; Results are from a state vector simulation
         n_shots (int): Number of times the circuit is run if ``from_samples=True``. Default: ``1000``
+        group_pauli_terms: Whether or not to group Pauli X/Y terms in the Hamiltonian together to reduce the measurement cost.
+            Default: ``None``; each of the Hamiltonian terms containing X/Y are in their own individual groups.
+        n_shots_per_pauli_term (bool): Whether or not ``n_shots`` is used for each Pauli term in the Hamiltonian, or for
+            *all* the terms in the Hamiltonian. Default: ``True``; ``n_shots`` are used to get the expectation value for each
+            term in the Hamiltonian.
+        shot_allocation: Iterable containing the number of shots to be allocated to each term in the Hamiltonian respectively if
+            n_shots_per_pauli_term is ``False``. Default: ``None``; shots are allocated based on the magnitudes of the coefficients
+            of the Hamiltonian terms.
 
     Returns:
-        Hamiltonian expectation value (float)
+        float: Hamiltonian expectation value
     """
-    if from_samples:
-        from functools import reduce
-
-        total = 0.0
-        # Iterate over each term in the Hamiltonian
-        for term in hamiltonian.terms:
-            # Get the target qubits and basis rotation gates from the Hamiltonian term
-            qubits = [int(factor.target_qubit) for factor in term.factors]
-            basis = [type(factor.gate) for factor in term.factors]
-            # Run a copy of the original circuit to get the output frequencies
+    if not from_samples:
+        # Expectation value from state vector simulation
+        result = circuit(nshots=1)
+        state_ket = result.state()
+        return hamiltonian.expectation(state_ket)
+
+    # From sample measurements:
+    # (Eventually) measurement_basis_rotations will be used to group up some terms so that one
+    # set of measurements can be used for multiple X/Y terms
+    grouped_terms = measurement_basis_rotations(hamiltonian, circuit.nqubits, grouping=group_pauli_terms)
+
+    # Check shot_allocation argument if not using n_shots_per_pauli_term
+    if not n_shots_per_pauli_term:
+        if shot_allocation is None:
+            shot_allocation = allocate_shots(grouped_terms, n_shots)
+        assert len(shot_allocation) == len(
+            grouped_terms
+        ), "shot_allocation list must be the same size as the number of grouped terms!"
+
+    total = 0.0
+    for _i, (measurement_gates, terms) in enumerate(grouped_terms):
+        if measurement_gates and terms:
             _circuit = circuit.copy()
-            _circuit.add(gates.M(*qubits, basis=basis))
-            result = _circuit(nshots=n_shots)
+            _circuit.add(measurement_gates)
+
+            # Number of shots used to run the circuit depends on n_shots_per_pauli_term
+            result = _circuit(nshots=n_shots if n_shots_per_pauli_term else shot_allocation[_i])
+
             frequencies = result.frequencies(binary=True)
-            # Only works for Z terms, raises an error if ham_term has X/Y terms
-            # total += SymbolicHamiltonian(
-            #              reduce(lambda x, y: x*y, term.factors, 1)
-            #          ).expectation_from_samples(frequencies, qubit_map=qubits)
-            # Workaround code to handle X/Y terms in the Hamiltonian:
-            # Get each Pauli string e.g. X0Y1
-            pauli = [factor.name for factor in term.factors]
-            # Replace each X and Y symbol with Z; then include the term coefficient
-            pauli_z = [Z(int(element[1:])) for element in pauli]
-            z_only_ham = SymbolicHamiltonian(term.coefficient * reduce(lambda x, y: x * y, pauli_z, 1.0))
-            # Can now apply expectation_from_samples directly
-            total += z_only_ham.expectation_from_samples(frequencies, qubit_map=qubits)
-        # Add the constant term if present. Note: Energies (in chemistry) are all real values
-        total += hamiltonian.constant.real
-        return total
-
-    # Expectation value from state vector simulation
-    result = circuit(nshots=1)
-    state_ket = result.state()
-    return hamiltonian.expectation(state_ket)
+            if frequencies:  # Needed because might have cases whereby no shots allocated to a group
+                # First term is all Z terms, can use expectation_from_samples directly.
+                # Otherwise, need to use the general pauli_term_measurement_expectation function
+                if _i > 0:
+                    total += sum(pauli_term_measurement_expectation(term, frequencies) for term in terms)
+                else:
+                    z_ham = SymbolicHamiltonian(sum(symbolic_term_to_symbol(term) for term in terms))
+                    total += z_ham.expectation_from_samples(frequencies)
+    # Add the constant term if present. Note: Energies (in chemistry) are all real values
+    total += hamiltonian.constant.real
+    return total

src/qibochem/measurement/optimization.py

@@ -0,0 +1,135 @@
+import numpy as np
+from qibo import gates
+
+
+def measurement_basis_rotations(hamiltonian, n_qubits, grouping=None):
+    """
+    Split up and sort the Hamiltonian terms to get the basis rotation gates to be applied to a quantum circuit for the
+    respective (group of) terms in the Hamiltonian
+
+    Args:
+        hamiltonian (SymbolicHamiltonian): Hamiltonian (that only contains X/Y terms?)
+        n_qubits: Number of qubits in the quantum circuit.
+        grouping: Whether or not to group the X/Y terms together, i.e. use the same set of measurements to get the expectation
+            values of a group of terms simultaneously. Default value of ``None`` will not group any terms together, which is
+            the only option currently implemented.
+
+    Returns:
+        list: List of two-tuples, with each tuple given as ([`list of measurement gates`], [term1, term2, ...]), where
+        term1, term2, ... are SymbolicTerms. The first tuple always corresponds to all the Z terms present, which will be two
+        empty lists - ``([], [])`` - if there are no Z terms present.
+    """
+    result = []
+    # Split up the Z and X/Y terms first
+    z_only_terms = [
+        term for term in hamiltonian.terms if not any(factor.name[0] in ("X", "Y") for factor in term.factors)
+    ]
+    xy_terms = [term for term in hamiltonian.terms if term not in z_only_terms]
+    # Add the Z terms into result first
+    if z_only_terms:
+        result.append(([gates.M(_i) for _i in range(n_qubits)], z_only_terms))
+    else:
+        result.append(([], []))
+    # Then add the X/Y terms in
+    if xy_terms:
+        if grouping is None:
+            result += [
+                (
+                    [
+                        gates.M(int(factor.target_qubit), basis=type(factor.gate))
+                        for factor in term.factors
+                        if factor.name[0] != "I"
+                    ],
+                    [
+                        term,
+                    ],
+                )
+                for term in xy_terms
+            ]
+        else:
+            raise NotImplementedError("Not ready yet!")
+    return result
+
+
+def allocate_shots(grouped_terms, n_shots, method=None, max_shots_per_term=None):
+    """
+    Allocate shots to each group of terms in the Hamiltonian for calculating the expectation value of the Hamiltonian.
+
+    Args:
+        grouped_terms (list): Output of measurement_basis_rotations(hamiltonian, n_qubits, grouping=None
+        n_shots (int): Total number of shots to be allocated
+        method (str): How to allocate the shots. The available options are: ``"c"``/``"coefficients"``: ``n_shots`` is distributed
+            based on the relative magnitudes of the term coefficients, ``"u"``/``"uniform"``: ``n_shots`` is distributed evenly
+            amongst each term. Default value: ``"c"``.
+        max_shots_per_term (int): Upper limit for the number of shots allocated to an individual group of terms. If not given,
+            will be defined as a fraction (largest coefficient over the sum of all coefficients in the Hamiltonian) of ``n_shots``.
+
+    Returns:
+        list: A list containing the number of shots to be used for each group of Pauli terms respectively.
+    """
+    if method is None:
+        method = "c"
+    if max_shots_per_term is None:
+        # Define based on the fraction of the term group with the largest coefficients w.r.t. sum of all coefficients.
+        term_coefficients = np.array(
+            [sum(abs(term.coefficient.real) for term in terms) for (_, terms) in grouped_terms]
+        )
+        max_shots_per_term = int(np.ceil(n_shots * (np.max(term_coefficients) / sum(term_coefficients))))
+        # max_shots_per_term = min(max_shots_per_term, 250)  #  Is there an optimal value - Explore further?
+    max_shots_per_term = min(n_shots, max_shots_per_term)  # Don't let max_shots_per_term > n_shots if manually defined
+
+    n_terms = len(grouped_terms)
+    shot_allocation = np.zeros(n_terms, dtype=int)
+
+    while True:
+        remaining_shots = n_shots - sum(shot_allocation)
+        if not remaining_shots:
+            break
+
+        # In case n_shots is too high s.t. max_shots_per_term is too low
+        # Increase max_shots_per_term, and redo the shot allocation
+        if np.min(shot_allocation) == max_shots_per_term:
+            max_shots_per_term = min(2 * max_shots_per_term, n_shots)
+            shot_allocation = np.zeros(n_terms, dtype=int)
+            continue
+
+        if method in ("c", "coefficients"):
+            # Split shots based on the relative magnitudes of the coefficients of the (group of) Pauli term(s)
+            # and only for terms that haven't reached the upper limit yet
+            term_coefficients = np.array(
+                [
+                    sum(abs(term.coefficient.real) for term in terms) if shots < max_shots_per_term else 0.0
+                    for shots, (_, terms) in zip(shot_allocation, grouped_terms)
+                ]
+            )
+
+            # Normalise term_coefficients, then get an initial distribution of remaining_shots
+            term_coefficients /= sum(term_coefficients)
+            _shot_allocation = (remaining_shots * term_coefficients).astype(int)
+            # Only keep the terms with >0 shots allocated, renormalise term_coefficients, and distribute again
+            term_coefficients *= _shot_allocation > 0
+            if _shot_allocation.any():
+                term_coefficients /= sum(term_coefficients)
+                _shot_allocation = (remaining_shots * term_coefficients).astype(int)
+            else:
+                # For distributing the remaining few shots, i.e. remaining_shots << n_terms
+                _shot_allocation = np.array(
+                    allocate_shots(grouped_terms, remaining_shots, max_shots_per_term=remaining_shots, method="u")
+                )
+
+        elif method in ("u", "uniform"):
+            # Uniform distribution of shots for every term. Extra shots are randomly distributed
+            _shot_allocation = np.array([remaining_shots // n_terms for _ in range(n_terms)])
+            if not _shot_allocation.any():
+                _shot_allocation = np.zeros(n_terms)
+                _shot_allocation[:remaining_shots] = 1
+                np.random.shuffle(_shot_allocation)
+
+        else:
+            raise NameError("Unknown method!")
+
+        # Add on to the current allocation, and set upper limit to the number of shots for a given term
+        shot_allocation += _shot_allocation.astype(int)
+        shot_allocation = np.clip(shot_allocation, 0, max_shots_per_term)
+
+    return shot_allocation.tolist()