|
| 1 | +# This code is part of a Qiskit project. |
| 2 | +# |
| 3 | +# (C) Copyright IBM 2022, 2024. |
| 4 | +# |
| 5 | +# This code is licensed under the Apache License, Version 2.0. You may |
| 6 | +# obtain a copy of this license in the LICENSE.txt file in the root directory |
| 7 | +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. |
| 8 | +# |
| 9 | +# Any modifications or derivative works of this code must retain this |
| 10 | +# copyright notice, and modified files need to carry a notice indicating |
| 11 | +# that they have been altered from the originals. |
| 12 | + |
| 13 | +"""Expectation value for a diagonal observable using a sampler primitive.""" |
| 14 | + |
| 15 | +from __future__ import annotations |
| 16 | + |
| 17 | +from collections.abc import Callable, Iterable, Mapping, MappingView, Sequence |
| 18 | +from typing import Any |
| 19 | + |
| 20 | +import numpy as np |
| 21 | +from qiskit.circuit import QuantumCircuit |
| 22 | +from qiskit.primitives import BaseEstimator, BaseSamplerV1, BaseSamplerV2 |
| 23 | +from qiskit.primitives.utils import _circuit_key, init_observable |
| 24 | +from qiskit.quantum_info import SparsePauliOp |
| 25 | +from qiskit.quantum_info.operators.base_operator import BaseOperator |
| 26 | +from qiskit_algorithms.algorithm_job import AlgorithmJob |
| 27 | +from qiskit_algorithms.minimum_eigensolvers.diagonal_estimator import _DiagonalEstimatorResult |
| 28 | + |
| 29 | + |
| 30 | +class _DiagonalEstimator(BaseEstimator): |
| 31 | + """An estimator for diagonal observables.""" |
| 32 | + |
| 33 | + def __init__( |
| 34 | + self, |
| 35 | + sampler: BaseSamplerV1 | BaseSamplerV2, |
| 36 | + aggregation: float | Callable[[Iterable[tuple[float, float]]], float] | None = None, |
| 37 | + callback: Callable[[Sequence[Mapping[str, Any]]], None] | None = None, |
| 38 | + **options, |
| 39 | + ) -> None: |
| 40 | + r"""Evaluate the expectation of quantum state with respect to a diagonal operator. |
| 41 | +
|
| 42 | + Args: |
| 43 | + sampler: The sampler used to evaluate the circuits. |
| 44 | + aggregation: The aggregation function to aggregate the measurement outcomes. If a float |
| 45 | + this specified the CVaR :math:`\alpha` parameter. |
| 46 | + callback: A callback which is given the best measurements of all circuits in each |
| 47 | + evaluation. |
| 48 | + run_options: Options for the sampler. |
| 49 | +
|
| 50 | + """ |
| 51 | + super().__init__(options=options) |
| 52 | + self._circuits: list[QuantumCircuit] = [] # See Qiskit pull request 11051 |
| 53 | + self._parameters: list[MappingView] = [] |
| 54 | + self._observables: list[SparsePauliOp] = [] |
| 55 | + |
| 56 | + self.sampler = sampler |
| 57 | + if not callable(aggregation): |
| 58 | + aggregation = _get_cvar_aggregation(aggregation) |
| 59 | + |
| 60 | + self.aggregation = aggregation |
| 61 | + self.callback = callback |
| 62 | + self._circuit_ids: dict[int, QuantumCircuit] = {} |
| 63 | + self._observable_ids: dict[int, BaseOperator] = {} |
| 64 | + |
| 65 | + def _run( |
| 66 | + self, |
| 67 | + circuits: Sequence[QuantumCircuit], |
| 68 | + observables: Sequence[BaseOperator], |
| 69 | + parameter_values: Sequence[Sequence[float]], |
| 70 | + **run_options, |
| 71 | + ) -> AlgorithmJob: |
| 72 | + circuit_indices = [] |
| 73 | + for circuit in circuits: |
| 74 | + key = _circuit_key(circuit) |
| 75 | + index = self._circuit_ids.get(key) |
| 76 | + if index is not None: |
| 77 | + circuit_indices.append(index) |
| 78 | + else: |
| 79 | + circuit_indices.append(len(self._circuits)) |
| 80 | + self._circuit_ids[key] = len(self._circuits) |
| 81 | + self._circuits.append(circuit) |
| 82 | + self._parameters.append(circuit.parameters) |
| 83 | + observable_indices = [] |
| 84 | + for observable in observables: |
| 85 | + index = self._observable_ids.get(id(observable)) |
| 86 | + if index is not None: |
| 87 | + observable_indices.append(index) |
| 88 | + else: |
| 89 | + observable_indices.append(len(self._observables)) |
| 90 | + self._observable_ids[id(observable)] = len(self._observables) |
| 91 | + converted_observable = init_observable(observable) |
| 92 | + _check_observable_is_diagonal(converted_observable) # check it's diagonal |
| 93 | + self._observables.append(converted_observable) |
| 94 | + job = AlgorithmJob( |
| 95 | + self._call, circuit_indices, observable_indices, parameter_values, **run_options |
| 96 | + ) |
| 97 | + job.submit() |
| 98 | + return job |
| 99 | + |
| 100 | + def _call( |
| 101 | + self, |
| 102 | + circuits: Sequence[int], |
| 103 | + observables: Sequence[int], |
| 104 | + parameter_values: Sequence[Sequence[float]], |
| 105 | + **run_options, |
| 106 | + ) -> _DiagonalEstimatorResult: |
| 107 | + if isinstance(self.sampler, BaseSamplerV1): |
| 108 | + job = self.sampler.run( |
| 109 | + [self._circuits[i] for i in circuits], |
| 110 | + parameter_values, |
| 111 | + **run_options, |
| 112 | + ) |
| 113 | + sampler_result = job.result() |
| 114 | + metadata = sampler_result.metadata |
| 115 | + samples = sampler_result.quasi_dists |
| 116 | + else: # BaseSamplerV2 |
| 117 | + job = self.sampler.run( |
| 118 | + [(self._circuits[i], val) for i, val in zip(circuits, parameter_values)], |
| 119 | + **run_options, |
| 120 | + ) |
| 121 | + sampler_pub_result = job.result() |
| 122 | + metadata = [] |
| 123 | + samples = [] |
| 124 | + for i, result in zip(circuits, sampler_pub_result): |
| 125 | + creg = self._circuits[i].cregs[0].name |
| 126 | + counts = getattr(result.data, creg).get_int_counts() |
| 127 | + shots = sum(counts.values()) |
| 128 | + samples.append({key: val / shots for key, val in counts.items()}) |
| 129 | + metadata.append(result.metadata) |
| 130 | + |
| 131 | + # a list of dictionaries containing: {state: (measurement probability, value)} |
| 132 | + evaluations: list[dict[int, tuple[float, float]]] = [ |
| 133 | + { |
| 134 | + state: (probability, _evaluate_sparsepauli(state, self._observables[i])) |
| 135 | + for state, probability in sampled.items() |
| 136 | + } |
| 137 | + for i, sampled in zip(observables, samples) |
| 138 | + ] |
| 139 | + |
| 140 | + results = np.array([self.aggregation(evaluated.values()) for evaluated in evaluations]) |
| 141 | + |
| 142 | + # get the best measurements |
| 143 | + best_measurements = [] |
| 144 | + num_qubits = self._circuits[0].num_qubits |
| 145 | + for evaluated in evaluations: |
| 146 | + best_result = min(evaluated.items(), key=lambda x: x[1][1]) |
| 147 | + best_measurements.append( |
| 148 | + { |
| 149 | + "state": best_result[0], |
| 150 | + "bitstring": bin(best_result[0])[2:].zfill(num_qubits), |
| 151 | + "value": best_result[1][1], |
| 152 | + "probability": best_result[1][0], |
| 153 | + } |
| 154 | + ) |
| 155 | + |
| 156 | + if self.callback is not None: |
| 157 | + self.callback(best_measurements) |
| 158 | + |
| 159 | + return _DiagonalEstimatorResult( |
| 160 | + values=results, metadata=metadata, best_measurements=best_measurements |
| 161 | + ) |
| 162 | + |
| 163 | + |
| 164 | +def _get_cvar_aggregation(alpha: float | None) -> Callable[[Iterable[tuple[float, float]]], float]: |
| 165 | + """Get the aggregation function for CVaR with confidence level ``alpha``.""" |
| 166 | + if alpha is None: |
| 167 | + alpha = 1 |
| 168 | + elif not 0 <= alpha <= 1: |
| 169 | + raise ValueError(f"alpha must be in [0, 1] but was {alpha}") |
| 170 | + |
| 171 | + # if alpha is close to 1 we can avoid the sorting |
| 172 | + if np.isclose(alpha, 1): |
| 173 | + |
| 174 | + def aggregate(measurements: Iterable[tuple[float, float]]) -> float: |
| 175 | + return sum(probability * value for probability, value in measurements) |
| 176 | + |
| 177 | + else: |
| 178 | + |
| 179 | + def aggregate(measurements: Iterable[tuple[float, float]]) -> float: |
| 180 | + # sort by values |
| 181 | + sorted_measurements = sorted(measurements, key=lambda x: x[1]) |
| 182 | + |
| 183 | + accumulated_percent = 0.0 # once alpha is reached, stop |
| 184 | + cvar = 0.0 |
| 185 | + for probability, value in sorted_measurements: |
| 186 | + cvar += value * min(probability, alpha - accumulated_percent) |
| 187 | + accumulated_percent += probability |
| 188 | + if accumulated_percent >= alpha: |
| 189 | + break |
| 190 | + |
| 191 | + return cvar / alpha |
| 192 | + |
| 193 | + return aggregate |
| 194 | + |
| 195 | + |
| 196 | +_PARITY = np.array([-1 if bin(i).count("1") % 2 else 1 for i in range(256)], dtype=np.complex128) |
| 197 | + |
| 198 | + |
| 199 | +def _evaluate_sparsepauli(state: int, observable: SparsePauliOp) -> float: |
| 200 | + packed_uint8 = np.packbits(observable.paulis.z, axis=1, bitorder="little") |
| 201 | + state_bytes = np.frombuffer(state.to_bytes(packed_uint8.shape[1], "little"), dtype=np.uint8) |
| 202 | + reduced = np.bitwise_xor.reduce(packed_uint8 & state_bytes, axis=1) |
| 203 | + return np.sum(observable.coeffs * _PARITY[reduced]) |
| 204 | + |
| 205 | + |
| 206 | +def _check_observable_is_diagonal(observable: SparsePauliOp) -> None: |
| 207 | + is_diagonal = not np.any(observable.paulis.x) |
| 208 | + if not is_diagonal: |
| 209 | + raise ValueError("The observable must be diagonal.") |
0 commit comments