From 940b94d6e6d0ce80df0b221a1587d1c4f6b8e090 Mon Sep 17 00:00:00 2001 From: Will Shanks Date: Thu, 14 Nov 2024 19:23:50 -0500 Subject: [PATCH 1/2] Adopt `qiskit.result.mitigation` into `qiskit_experiments.data_processing` (#1484) The readout error mitigator classes `LocalReadoutMitigator` and `CorrelatedReadoutMitigator` have been copied into `qiskit_experiments.data_processing.mitigation` from `qiskit.result.mitigation` in order to prepare for the deprecation of these classes from Qiskit. The experiments that previously used the Qiskit versions of these mitigator classes now use the Experiments versions. The manual has also been updated to use the new class locations. Additionally, the readout mitigation manual was refactored to avoid using private functions and methods. Some adjacent `jupyter-execute` cells were merged because in the rendered docs they look like one cell any way, but they still have individual "copy to clipboard" buttons which is confusing. You think clicking the button will get the merged cell contents but it only gets a section of it. --- .../measurement/readout_mitigation.rst | 47 +- .../data_processing/__init__.py | 12 + .../data_processing/mitigation/__init__.py | 22 + .../mitigation/base_readout_mitigator.py | 81 +++ .../correlated_readout_mitigator.py | 270 ++++++++++ .../mitigation/local_readout_mitigator.py | 320 ++++++++++++ .../data_processing/mitigation/utils.py | 163 ++++++ .../correlated_readout_error_analysis.py | 2 +- .../analysis/local_readout_error_analysis.py | 2 +- .../library/tomography/basis/pauli_basis.py | 2 +- .../adopt-mitigators-d7ad0e2f3cd2fa57.yaml | 18 + test/data_processing/test_mitigators.py | 477 ++++++++++++++++++ .../tomography/test_process_tomography.py | 2 +- .../tomography/test_state_tomography.py | 2 +- 14 files changed, 1388 insertions(+), 32 deletions(-) create mode 100644 qiskit_experiments/data_processing/mitigation/__init__.py create mode 100644 qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py create mode 100644 qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py create mode 100644 qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py create mode 100644 qiskit_experiments/data_processing/mitigation/utils.py create mode 100644 releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml create mode 100644 test/data_processing/test_mitigators.py diff --git a/docs/manuals/measurement/readout_mitigation.rst b/docs/manuals/measurement/readout_mitigation.rst index c0672a6118..80e74f65fa 100644 --- a/docs/manuals/measurement/readout_mitigation.rst +++ b/docs/manuals/measurement/readout_mitigation.rst @@ -47,25 +47,16 @@ experiments to generate the corresponding mitigators. import numpy as np import matplotlib.pyplot as plt from qiskit import QuantumCircuit + from qiskit.quantum_info import Operator from qiskit.visualization import plot_distribution + from qiskit_experiments.data_processing import LocalReadoutMitigator from qiskit_experiments.library import LocalReadoutError, CorrelatedReadoutError from qiskit_aer import AerSimulator from qiskit_ibm_runtime.fake_provider import FakePerth - from qiskit.result.mitigation.utils import ( - expval_with_stddev, - str2diag, - counts_probability_vector - ) - backend = AerSimulator.from_backend(FakePerth()) -.. jupyter-execute:: - - shots = 1024 - qubits = [0,1,2,3] - num_qubits = len(qubits) Standard mitigation experiment ------------------------------ @@ -76,13 +67,14 @@ circuits, one for all “0” and one for all “1” results. .. jupyter-execute:: + shots = 1024 + qubits = [0,1,2,3] + num_qubits = len(qubits) + exp = LocalReadoutError(qubits) for c in exp.circuits(): print(c) - -.. jupyter-execute:: - exp.analysis.set_options(plot=True) result = exp.run(backend) mitigator = result.analysis_results("Local Readout Mitigator").value @@ -102,9 +94,9 @@ The individual mitigation matrices can be read off the mitigator. .. jupyter-execute:: - for m in mitigator._mitigation_mats: - print(m) - print() + for qubit in mitigator.qubits: + print(f"Qubit: {qubit}") + print(mitigator.mitigation_matrix(qubits=qubit)) Mitigation example @@ -118,13 +110,9 @@ Mitigation example qc.cx(i - 1, i) qc.measure_all() -.. jupyter-execute:: - counts = backend.run(qc, shots=shots, seed_simulator=42, method="density_matrix").result().get_counts() unmitigated_probs = {label: count / shots for label, count in counts.items()} -.. jupyter-execute:: - mitigated_quasi_probs = mitigator.quasi_probabilities(counts) mitigated_stddev = mitigated_quasi_probs._stddev_upper_bound mitigated_probs = (mitigated_quasi_probs.nearest_probability_distribution().binary_probabilities()) @@ -144,15 +132,20 @@ Expectation value .. jupyter-execute:: diagonal_labels = ["ZZZZ", "ZIZI", "IZII", "1ZZ0"] - ideal_expectation = [] - diagonals = [str2diag(d) for d in diagonal_labels] + diagonals = [ + np.diag(np.real(Operator.from_label(d).to_matrix())) + for d in diagonal_labels + ] + + # Create a mitigator with no mitigation so that we can use its + # expectation_values method to generate an unmitigated expectation value to + # compare to the mitigated one. + identity_mitigator = LocalReadoutMitigator([np.eye(2) for _ in range(4)]) + qubit_index = {i: i for i in range(num_qubits)} - unmitigated_probs_vector, _ = counts_probability_vector(unmitigated_probs, qubit_index=qubit_index) - unmitigated_expectation = [expval_with_stddev(d, unmitigated_probs_vector, shots) for d in diagonals] + unmitigated_expectation = [identity_mitigator.expectation_value(counts, d) for d in diagonals] mitigated_expectation = [mitigator.expectation_value(counts, d) for d in diagonals] -.. jupyter-execute:: - mitigated_expectation_values, mitigated_stddev = zip(*mitigated_expectation) unmitigated_expectation_values, unmitigated_stddev = zip(*unmitigated_expectation) legend = ['Mitigated Expectation', 'Unmitigated Expectation'] diff --git a/qiskit_experiments/data_processing/__init__.py b/qiskit_experiments/data_processing/__init__.py index 56a05ddfa3..b075c1effe 100644 --- a/qiskit_experiments/data_processing/__init__.py +++ b/qiskit_experiments/data_processing/__init__.py @@ -82,6 +82,15 @@ BaseDiscriminator SkLDA SkQDA + +Mitigators +========== +.. autosummary:: + :toctree: ../stubs/ + + BaseReadoutMitigator + LocalReadoutMitigator + CorrelatedReadoutMitigator """ from .data_action import DataAction, TrainableDataAction @@ -104,4 +113,7 @@ from .data_processor import DataProcessor from .discriminator import BaseDiscriminator +from .mitigation.base_readout_mitigator import BaseReadoutMitigator +from .mitigation.correlated_readout_mitigator import CorrelatedReadoutMitigator +from .mitigation.local_readout_mitigator import LocalReadoutMitigator from .sklearn_discriminators import SkLDA, SkQDA diff --git a/qiskit_experiments/data_processing/mitigation/__init__.py b/qiskit_experiments/data_processing/mitigation/__init__.py new file mode 100644 index 0000000000..5cdc4ebf1d --- /dev/null +++ b/qiskit_experiments/data_processing/mitigation/__init__.py @@ -0,0 +1,22 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2021. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Readout error mitigation.""" +from .base_readout_mitigator import BaseReadoutMitigator +from .correlated_readout_mitigator import CorrelatedReadoutMitigator +from .local_readout_mitigator import LocalReadoutMitigator +from .utils import ( + counts_probability_vector, + expval_with_stddev, + stddev, + str2diag, +) diff --git a/qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py b/qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py new file mode 100644 index 0000000000..c37fa333ea --- /dev/null +++ b/qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py @@ -0,0 +1,81 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Base class for readout error mitigation. +""" + +from abc import ABC, abstractmethod +from typing import Optional, List, Iterable, Tuple, Union, Callable + +import numpy as np + +from qiskit.result.counts import Counts +from qiskit.result.distributions.quasi import QuasiDistribution + + +class BaseReadoutMitigator(ABC): + """Base readout error mitigator class.""" + + @abstractmethod + def quasi_probabilities( + self, + data: Counts, + qubits: Iterable[int] = None, + clbits: Optional[List[int]] = None, + shots: Optional[int] = None, + ) -> QuasiDistribution: + """Convert counts to a dictionary of quasi-probabilities + + Args: + data: Counts to be mitigated. + qubits: the physical qubits measured to obtain the counts clbits. + If None these are assumed to be qubits [0, ..., N-1] + for N-bit counts. + clbits: Optional, marginalize counts to just these bits. + shots: Optional, the total number of shots, if None shots will + be calculated as the sum of all counts. + + Returns: + QuasiDistribution: A dictionary containing pairs of [output, mean] where "output" + is the key in the dictionaries, + which is the length-N bitstring of a measured standard basis state, + and "mean" is the mean of non-zero quasi-probability estimates. + """ + + @abstractmethod + def expectation_value( + self, + data: Counts, + diagonal: Union[Callable, dict, str, np.ndarray], + qubits: Iterable[int] = None, + clbits: Optional[List[int]] = None, + shots: Optional[int] = None, + ) -> Tuple[float, float]: + """Calculate the expectation value of a diagonal Hermitian operator. + + Args: + data: Counts object to be mitigated. + diagonal: the diagonal operator. This may either be specified + as a string containing I,Z,0,1 characters, or as a + real valued 1D array_like object supplying the full diagonal, + or as a dictionary, or as Callable. + qubits: the physical qubits measured to obtain the counts clbits. + If None these are assumed to be qubits [0, ..., N-1] + for N-bit counts. + clbits: Optional, marginalize counts to just these bits. + shots: Optional, the total number of shots, if None shots will + be calculated as the sum of all counts. + + Returns: + The mean and an upper bound of the standard deviation of operator + expectation value calculated from the current counts. + """ diff --git a/qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py b/qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py new file mode 100644 index 0000000000..8099347285 --- /dev/null +++ b/qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py @@ -0,0 +1,270 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Readout mitigator class based on the A-matrix inversion method +""" + +import math +from typing import Optional, List, Tuple, Iterable, Callable, Union, Dict +import numpy as np + +from qiskit.exceptions import QiskitError + +from qiskit.result.distributions.quasi import QuasiDistribution +from qiskit.result.counts import Counts +from .base_readout_mitigator import BaseReadoutMitigator +from .utils import counts_probability_vector, z_diagonal, str2diag + + +class CorrelatedReadoutMitigator(BaseReadoutMitigator): + """N-qubit readout error mitigator. + + Mitigates :meth:`expectation_value` and :meth:`quasi_probabilities`. + The mitigation_matrix should be calibrated using qiskit experiments. + This mitigation method should be used in case the readout errors of the qubits + are assumed to be correlated. The mitigation_matrix of *N* qubits is of size + :math:`2^N x 2^N` so the mitigation complexity is :math:`O(4^N)`. + """ + + def __init__(self, assignment_matrix: np.ndarray, qubits: Optional[Iterable[int]] = None): + """Initialize a CorrelatedReadoutMitigator + + Args: + assignment_matrix: readout error assignment matrix. + qubits: Optional, the measured physical qubits for mitigation. + + Raises: + QiskitError: matrix size does not agree with number of qubits + """ + if np.any(assignment_matrix < 0) or not np.allclose(np.sum(assignment_matrix, axis=0), 1): + raise QiskitError("Assignment matrix columns must be valid probability distributions") + assignment_matrix = np.asarray(assignment_matrix, dtype=float) + matrix_qubits_num = int(math.log2(assignment_matrix.shape[0])) + if qubits is None: + self._num_qubits = matrix_qubits_num + self._qubits = range(self._num_qubits) + else: + if len(qubits) != matrix_qubits_num: + raise QiskitError( + f"The number of given qubits ({len(qubits)}) is different than the number of " + f"qubits inferred from the matrices ({matrix_qubits_num})" + ) + self._qubits = qubits + self._num_qubits = len(self._qubits) + self._qubit_index = dict(zip(self._qubits, range(self._num_qubits))) + self._assignment_mat = assignment_matrix + self._mitigation_mats = {} + + @property + def settings(self) -> Dict: + """Return settings.""" + return {"assignment_matrix": self._assignment_mat, "qubits": self._qubits} + + def expectation_value( + self, + data: Counts, + diagonal: Union[Callable, dict, str, np.ndarray] = None, + qubits: Iterable[int] = None, + clbits: Optional[List[int]] = None, + shots: Optional[int] = None, + ) -> Tuple[float, float]: + r"""Compute the mitigated expectation value of a diagonal observable. + + This computes the mitigated estimator of + :math:`\langle O \rangle = \mbox{Tr}[\rho. O]` of a diagonal observable + :math:`O = \sum_{x\in\{0, 1\}^n} O(x)|x\rangle\!\langle x|`. + + Args: + data: Counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the default value is + :math:`[1, -1]^\otimes n`. + qubits: Optional, the measured physical qubits the count + bitstrings correspond to. If None qubits are assumed to be + :math:`[0, ..., n-1]`. + clbits: Optional, if not None marginalize counts to the specified bits. + shots: the number of shots. + + Returns: + (float, float): the expectation value and an upper bound of the standard deviation. + + Additional Information: + The diagonal observable :math:`O` is input using the ``diagonal`` kwarg as + a list or Numpy array :math:`[O(0), ..., O(2^n -1)]`. If no diagonal is specified + the diagonal of the Pauli operator + :math`O = \mbox{diag}(Z^{\otimes n}) = [1, -1]^{\otimes n}` is used. + The ``clbits`` kwarg is used to marginalize the input counts dictionary + over the specified bit-values, and the ``qubits`` kwarg is used to specify + which physical qubits these bit-values correspond to as + ``circuit.measure(qubits, clbits)``. + """ + + if qubits is None: + qubits = self._qubits + probs_vec, shots = counts_probability_vector( + data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits + ) + + # Get qubit mitigation matrix and mitigate probs + mit_mat = self.mitigation_matrix(qubits) + + # Get operator coeffs + if diagonal is None: + diagonal = z_diagonal(2**self._num_qubits) + elif isinstance(diagonal, str): + diagonal = str2diag(diagonal) + + # Apply transpose of mitigation matrix + coeffs = mit_mat.T.dot(diagonal) + expval = coeffs.dot(probs_vec) + stddev_upper_bound = self.stddev_upper_bound(shots) + + return (expval, stddev_upper_bound) + + def quasi_probabilities( + self, + data: Counts, + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, + shots: Optional[int] = None, + ) -> QuasiDistribution: + """Compute mitigated quasi probabilities value. + + Args: + data: counts object + qubits: qubits the count bitstrings correspond to. + clbits: Optional, marginalize counts to just these bits. + shots: Optional, the total number of shots, if None shots will + be calculated as the sum of all counts. + + Returns: + QuasiDistribution: A dictionary containing pairs of [output, mean] where "output" + is the key in the dictionaries, + which is the length-N bitstring of a measured standard basis state, + and "mean" is the mean of non-zero quasi-probability estimates. + """ + if qubits is None: + qubits = self._qubits + probs_vec, calculated_shots = counts_probability_vector( + data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits + ) + if shots is None: + shots = calculated_shots + + # Get qubit mitigation matrix and mitigate probs + mit_mat = self.mitigation_matrix(qubits) + + # Apply transpose of mitigation matrix + probs_vec = mit_mat.dot(probs_vec) + probs_dict = {} + for index, _ in enumerate(probs_vec): + probs_dict[index] = probs_vec[index] + + quasi_dist = QuasiDistribution( + probs_dict, stddev_upper_bound=self.stddev_upper_bound(shots) + ) + + return quasi_dist + + def mitigation_matrix(self, qubits: List[int] = None) -> np.ndarray: + r"""Return the readout mitigation matrix for the specified qubits. + + The mitigation matrix :math:`A^{-1}` is defined as the inverse of the + :meth:`assignment_matrix` :math:`A`. + + Args: + qubits: Optional, qubits being measured. + + Returns: + np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`. + """ + if qubits is None: + qubits = self._qubits + qubits = tuple(sorted(qubits)) + + # Check for cached mitigation matrix + # if not present compute + if qubits not in self._mitigation_mats: + marginal_matrix = self.assignment_matrix(qubits) + try: + mit_mat = np.linalg.inv(marginal_matrix) + except np.linalg.LinAlgError: + # Use pseudo-inverse if matrix is singular + mit_mat = np.linalg.pinv(marginal_matrix) + self._mitigation_mats[qubits] = mit_mat + + return self._mitigation_mats[qubits] + + def assignment_matrix(self, qubits: List[int] = None) -> np.ndarray: + r"""Return the readout assignment matrix for specified qubits. + + The assignment matrix is the stochastic matrix :math:`A` which assigns + a noisy readout probability distribution to an ideal input + readout distribution: :math:`P(i|j) = \langle i | A | j \rangle`. + + Args: + qubits: Optional, qubits being measured. + + Returns: + np.ndarray: the assignment matrix A. + """ + if qubits is None: + qubits = self._qubits + if qubits == self._num_qubits: + return self._assignment_mat + + if isinstance(qubits, int): + qubits = [qubits] + + qubit_indices = [self._qubit_index[qubit] for qubit in qubits] + # Compute marginal matrix + axis = tuple( + self._num_qubits - 1 - i for i in set(range(self._num_qubits)).difference(qubit_indices) + ) + num_qubits = len(qubits) + + new_amat = np.zeros(2 * [2**num_qubits], dtype=float) + for i, col in enumerate(self._assignment_mat.T[self._keep_indexes(qubit_indices)]): + new_amat[i] = ( + np.reshape(col, self._num_qubits * [2]).sum(axis=axis).reshape([2**num_qubits]) + ) + new_amat = new_amat.T + return new_amat + + @staticmethod + def _keep_indexes(qubits): + indexes = [0] + for i in sorted(qubits): + indexes += [idx + (1 << i) for idx in indexes] + return indexes + + def _compute_gamma(self): + """Compute gamma for N-qubit mitigation""" + mitmat = self.mitigation_matrix(qubits=self._qubits) + return np.max(np.sum(np.abs(mitmat), axis=0)) + + def stddev_upper_bound(self, shots: int): + """Return an upper bound on standard deviation of expval estimator. + + Args: + shots: Number of shots used for expectation value measurement. + + Returns: + float: the standard deviation upper bound. + """ + gamma = self._compute_gamma() + return gamma / math.sqrt(shots) + + @property + def qubits(self) -> Tuple[int]: + """The device qubits for this mitigator""" + return self._qubits diff --git a/qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py b/qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py new file mode 100644 index 0000000000..4e843dbbd4 --- /dev/null +++ b/qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py @@ -0,0 +1,320 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Readout mitigator class based on the 1-qubit local tensored mitigation method +""" + + +import math +from typing import Optional, List, Tuple, Iterable, Callable, Union, Dict +import numpy as np + +from qiskit.exceptions import QiskitError +from qiskit.result.distributions.quasi import QuasiDistribution +from qiskit.result.counts import Counts +from .base_readout_mitigator import BaseReadoutMitigator +from .utils import counts_probability_vector, z_diagonal, str2diag + + +class LocalReadoutMitigator(BaseReadoutMitigator): + """1-qubit tensor product readout error mitigator. + + Mitigates :meth:`expectation_value` and :meth:`quasi_probabilities`. + The mitigator should either be calibrated using qiskit experiments, + or calculated directly from the backend properties. + This mitigation method should be used in case the readout errors of the qubits + are assumed to be uncorrelated. For *N* qubits there are *N* mitigation matrices, + each of size :math:`2 x 2` and the mitigation complexity is :math:`O(2^N)`, + so it is more efficient than the :class:`CorrelatedReadoutMitigator` class. + """ + + def __init__( + self, + assignment_matrices: Optional[List[np.ndarray]] = None, + qubits: Optional[Iterable[int]] = None, + backend=None, + ): + """Initialize a LocalReadoutMitigator + + Args: + assignment_matrices: Optional, list of single-qubit readout error assignment matrices. + qubits: Optional, the measured physical qubits for mitigation. + backend: Optional, backend name. + + Raises: + QiskitError: matrices sizes do not agree with number of qubits + """ + if assignment_matrices is None: + assignment_matrices = self._from_backend(backend, qubits) + else: + assignment_matrices = [np.asarray(amat, dtype=float) for amat in assignment_matrices] + for amat in assignment_matrices: + if np.any(amat < 0) or not np.allclose(np.sum(amat, axis=0), 1): + raise QiskitError( + "Assignment matrix columns must be valid probability distributions" + ) + if qubits is None: + self._num_qubits = len(assignment_matrices) + self._qubits = range(self._num_qubits) + else: + if len(qubits) != len(assignment_matrices): + raise QiskitError( + f"The number of given qubits ({len(qubits)}) is different than the number of qubits " + f"inferred from the matrices ({len(assignment_matrices)})" + ) + self._qubits = qubits + self._num_qubits = len(self._qubits) + + self._qubit_index = dict(zip(self._qubits, range(self._num_qubits))) + self._assignment_mats = assignment_matrices + self._mitigation_mats = np.zeros([self._num_qubits, 2, 2], dtype=float) + self._gammas = np.zeros(self._num_qubits, dtype=float) + + for i in range(self._num_qubits): + mat = self._assignment_mats[i] + # Compute Gamma values + error0 = mat[1, 0] + error1 = mat[0, 1] + self._gammas[i] = (1 + abs(error0 - error1)) / (1 - error0 - error1) + # Compute inverse mitigation matrix + try: + ainv = np.linalg.inv(mat) + except np.linalg.LinAlgError: + ainv = np.linalg.pinv(mat) + self._mitigation_mats[i] = ainv + + @property + def settings(self) -> Dict: + """Return settings.""" + return {"assignment_matrices": self._assignment_mats, "qubits": self._qubits} + + def expectation_value( + self, + data: Counts, + diagonal: Union[Callable, dict, str, np.ndarray] = None, + qubits: Iterable[int] = None, + clbits: Optional[List[int]] = None, + shots: Optional[int] = None, + ) -> Tuple[float, float]: + r"""Compute the mitigated expectation value of a diagonal observable. + + This computes the mitigated estimator of + :math:`\langle O \rangle = \mbox{Tr}[\rho. O]` of a diagonal observable + :math:`O = \sum_{x\in\{0, 1\}^n} O(x)|x\rangle\!\langle x|`. + + Args: + data: Counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the default value is + :math:`[1, -1]^\otimes n`. + qubits: Optional, the measured physical qubits the count + bitstrings correspond to. If None qubits are assumed to be + :math:`[0, ..., n-1]`. + clbits: Optional, if not None marginalize counts to the specified bits. + shots: the number of shots. + + Returns: + (float, float): the expectation value and an upper bound of the standard deviation. + + Additional Information: + The diagonal observable :math:`O` is input using the ``diagonal`` kwarg as + a list or Numpy array :math:`[O(0), ..., O(2^n -1)]`. If no diagonal is specified + the diagonal of the Pauli operator + :math`O = \mbox{diag}(Z^{\otimes n}) = [1, -1]^{\otimes n}` is used. + The ``clbits`` kwarg is used to marginalize the input counts dictionary + over the specified bit-values, and the ``qubits`` kwarg is used to specify + which physical qubits these bit-values correspond to as + ``circuit.measure(qubits, clbits)``. + """ + if qubits is None: + qubits = self._qubits + num_qubits = len(qubits) + probs_vec, shots = counts_probability_vector( + data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits + ) + + # Get qubit mitigation matrix and mitigate probs + qubit_indices = [self._qubit_index[qubit] for qubit in qubits] + ainvs = self._mitigation_mats[qubit_indices] + + # Get operator coeffs + if diagonal is None: + diagonal = z_diagonal(2**num_qubits) + elif isinstance(diagonal, str): + diagonal = str2diag(diagonal) + + # Apply transpose of mitigation matrix + coeffs = np.reshape(diagonal, num_qubits * [2]) + einsum_args = [coeffs, list(range(num_qubits))] + for i, ainv in enumerate(reversed(ainvs)): + einsum_args += [ainv.T, [num_qubits + i, i]] + einsum_args += [list(range(num_qubits, 2 * num_qubits))] + coeffs = np.einsum(*einsum_args).ravel() + + expval = coeffs.dot(probs_vec) + stddev_upper_bound = self.stddev_upper_bound(shots, qubits) + + return (expval, stddev_upper_bound) + + def quasi_probabilities( + self, + data: Counts, + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, + shots: Optional[int] = None, + ) -> QuasiDistribution: + """Compute mitigated quasi probabilities value. + + Args: + data: counts object + qubits: qubits the count bitstrings correspond to. + clbits: Optional, marginalize counts to just these bits. + shots: Optional, the total number of shots, if None shots will + be calculated as the sum of all counts. + + Returns: + QuasiDistribution: A dictionary containing pairs of [output, mean] where "output" + is the key in the dictionaries, + which is the length-N bitstring of a measured standard basis state, + and "mean" is the mean of non-zero quasi-probability estimates. + + Raises: + QiskitError: if qubit and clbit kwargs are not valid. + """ + if qubits is None: + qubits = self._qubits + + num_qubits = len(qubits) + + probs_vec, calculated_shots = counts_probability_vector( + data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits + ) + if shots is None: + shots = calculated_shots + + # Get qubit mitigation matrix and mitigate probs + qubit_indices = [self._qubit_index[qubit] for qubit in qubits] + ainvs = self._mitigation_mats[qubit_indices] + + # Apply transpose of mitigation matrix + prob_tens = np.reshape(probs_vec, num_qubits * [2]) + einsum_args = [prob_tens, list(range(num_qubits))] + for i, ainv in enumerate(reversed(ainvs)): + einsum_args += [ainv, [num_qubits + i, i]] + einsum_args += [list(range(num_qubits, 2 * num_qubits))] + probs_vec = np.einsum(*einsum_args).ravel() + + probs_dict = {} + for index, _ in enumerate(probs_vec): + probs_dict[index] = probs_vec[index] + + quasi_dist = QuasiDistribution( + probs_dict, shots=shots, stddev_upper_bound=self.stddev_upper_bound(shots, qubits) + ) + return quasi_dist + + def mitigation_matrix(self, qubits: Optional[Union[List[int], int]] = None) -> np.ndarray: + r"""Return the measurement mitigation matrix for the specified qubits. + + The mitigation matrix :math:`A^{-1}` is defined as the inverse of the + :meth:`assignment_matrix` :math:`A`. + + Args: + qubits: Optional, qubits being measured for operator expval. + if a single int is given, it is assumed to be the index + of the qubit in self._qubits + + Returns: + np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`. + """ + if qubits is None: + qubits = self._qubits + if isinstance(qubits, int): + qubits = [self._qubits[qubits]] + qubit_indices = [self._qubit_index[qubit] for qubit in qubits] + mat = self._mitigation_mats[qubit_indices[0]] + for i in qubit_indices[1:]: + mat = np.kron(self._mitigation_mats[i], mat) + return mat + + def assignment_matrix(self, qubits: List[int] = None) -> np.ndarray: + r"""Return the measurement assignment matrix for specified qubits. + + The assignment matrix is the stochastic matrix :math:`A` which assigns + a noisy measurement probability distribution to an ideal input + measurement distribution: :math:`P(i|j) = \langle i | A | j \rangle`. + + Args: + qubits: Optional, qubits being measured for operator expval. + + Returns: + np.ndarray: the assignment matrix A. + """ + if qubits is None: + qubits = self._qubits + if isinstance(qubits, int): + qubits = [qubits] + qubit_indices = [self._qubit_index[qubit] for qubit in qubits] + mat = self._assignment_mats[qubit_indices[0]] + for i in qubit_indices[1:]: + mat = np.kron(self._assignment_mats[i], mat) + return mat + + def _compute_gamma(self, qubits=None): + """Compute gamma for N-qubit mitigation""" + if qubits is None: + gammas = self._gammas + else: + qubit_indices = [self._qubit_index[qubit] for qubit in qubits] + gammas = self._gammas[qubit_indices] + return np.prod(gammas) + + def stddev_upper_bound(self, shots: int, qubits: List[int] = None): + """Return an upper bound on standard deviation of expval estimator. + + Args: + shots: Number of shots used for expectation value measurement. + qubits: qubits being measured for operator expval. + + Returns: + float: the standard deviation upper bound. + """ + gamma = self._compute_gamma(qubits=qubits) + return gamma / math.sqrt(shots) + + def _from_backend(self, backend, qubits): + """Calculates amats from backend properties readout_error""" + backend_qubits = backend.properties().qubits + if qubits is not None: + if any(qubit >= len(backend_qubits) for qubit in qubits): + raise QiskitError("The chosen backend does not contain the specified qubits.") + reduced_backend_qubits = [backend_qubits[i] for i in qubits] + backend_qubits = reduced_backend_qubits + num_qubits = len(backend_qubits) + + amats = np.zeros([num_qubits, 2, 2], dtype=float) + + for qubit_idx, qubit_prop in enumerate(backend_qubits): + for prop in qubit_prop: + if prop.name == "prob_meas0_prep1": + (amats[qubit_idx])[0, 1] = prop.value + (amats[qubit_idx])[1, 1] = 1 - prop.value + if prop.name == "prob_meas1_prep0": + (amats[qubit_idx])[1, 0] = prop.value + (amats[qubit_idx])[0, 0] = 1 - prop.value + + return amats + + @property + def qubits(self) -> Tuple[int]: + """The device qubits for this mitigator""" + return self._qubits diff --git a/qiskit_experiments/data_processing/mitigation/utils.py b/qiskit_experiments/data_processing/mitigation/utils.py new file mode 100644 index 0000000000..6f4b9d7871 --- /dev/null +++ b/qiskit_experiments/data_processing/mitigation/utils.py @@ -0,0 +1,163 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Readout mitigation data handling utils + +.. note:: + + None of these functions are intended for public use. +""" + +import logging +import math +from typing import Optional, List, Tuple, Dict +import numpy as np + +from qiskit.exceptions import QiskitError +from qiskit.result import Counts, marginal_counts + +logger = logging.getLogger(__name__) + + +def z_diagonal(dim, dtype=float): + r"""Return the diagonal for the operator :math:`Z^\otimes n`""" + parity = np.zeros(dim, dtype=dtype) + for i in range(dim): + parity[i] = bin(i)[2:].count("1") + return (-1) ** np.mod(parity, 2) + + +def expval_with_stddev(coeffs: np.ndarray, probs: np.ndarray, shots: int) -> Tuple[float, float]: + """Compute expectation value and standard deviation. + Args: + coeffs: array of diagonal operator coefficients. + probs: array of measurement probabilities. + shots: total number of shots to obtain probabilities. + Returns: + tuple: (expval, stddev) expectation value and standard deviation. + """ + # Compute expval + expval = coeffs.dot(probs) + + # Compute variance + sq_expval = (coeffs**2).dot(probs) + variance = (sq_expval - expval**2) / shots + + # Compute standard deviation + if variance < 0 and not np.isclose(variance, 0): + logger.warning( + "Encountered a negative variance in expectation value calculation." + "(%f). Setting standard deviation of result to 0.", + variance, + ) + calc_stddev = math.sqrt(variance) if variance > 0 else 0.0 + return [expval, calc_stddev] + + +def stddev(probs, shots): + """Calculate stddev dict""" + ret = {} + for key, prob in probs.items(): + std_err = math.sqrt(prob * (1 - prob) / shots) + ret[key] = std_err + return ret + + +def str2diag(string): + """Transform diagonal from a string to a numpy array""" + chars = { + "I": np.array([1, 1], dtype=float), + "Z": np.array([1, -1], dtype=float), + "0": np.array([1, 0], dtype=float), + "1": np.array([0, 1], dtype=float), + } + ret = np.array([1], dtype=float) + for i in reversed(string): + if i not in chars: + raise QiskitError(f"Invalid diagonal string character {i}") + ret = np.kron(chars[i], ret) + return ret + + +def counts_to_vector(counts: Counts, num_qubits: int) -> Tuple[np.ndarray, int]: + """Transforms Counts to a probability vector""" + vec = np.zeros(2**num_qubits, dtype=float) + shots = 0 + for key, val in counts.items(): + shots += val + vec[int(key, 2)] = val + vec /= shots + return vec, shots + + +def remap_qubits( + vec: np.ndarray, num_qubits: int, qubits: Optional[List[int]] = None +) -> np.ndarray: + """Remapping the qubits""" + if qubits is not None: + if len(qubits) != num_qubits: + raise QiskitError("Num qubits does not match vector length.") + axes = [num_qubits - 1 - i for i in reversed(np.argsort(qubits))] + vec = np.reshape(vec, num_qubits * [2]).transpose(axes).reshape(vec.shape) + return vec + + +def marginalize_counts( + counts: Counts, + qubit_index: Dict[int, int], + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, +) -> np.ndarray: + """Marginalization of the Counts. Verify that number of clbits equals to the number of qubits.""" + if clbits is not None: + qubits_len = len(qubits) if not qubits is None else 0 + clbits_len = len(clbits) if not clbits is None else 0 + if clbits_len not in (0, qubits_len): + raise QiskitError( + f"Num qubits ({qubits_len}) does not match number of clbits ({clbits_len})." + ) + counts = marginal_counts(counts, clbits) + if clbits is None and qubits is not None: + clbits = [qubit_index[qubit] for qubit in qubits] + counts = marginal_counts(counts, clbits) + return counts + + +def counts_probability_vector( + counts: Counts, + qubit_index: Dict[int, int], + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, +) -> Tuple[np.ndarray, int]: + """Compute a probability vector for all count outcomes. + + Args: + counts: counts object + qubit_index: For each qubit, its index in the mitigator qubits list + qubits: qubits the count bitstrings correspond to. + clbits: Optional, marginalize counts to just these bits. + + Raises: + QiskitError: if qubits and clbits kwargs are not valid. + + Returns: + np.ndarray: a probability vector for all count outcomes. + int: Number of shots in the counts + """ + counts = marginalize_counts(counts, qubit_index, qubits, clbits) + if qubits is not None: + num_qubits = len(qubits) + else: + num_qubits = len(qubit_index.keys()) + vec, shots = counts_to_vector(counts, num_qubits) + vec = remap_qubits(vec, num_qubits, qubits) + return vec, shots diff --git a/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py b/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py index 73dece5171..796c375d9f 100644 --- a/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py @@ -15,7 +15,7 @@ from typing import List, Tuple import numpy as np import matplotlib.pyplot as plt -from qiskit.result import CorrelatedReadoutMitigator +from qiskit_experiments.data_processing import CorrelatedReadoutMitigator from qiskit_experiments.framework import ExperimentData from qiskit_experiments.framework.matplotlib import get_non_gui_ax from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options diff --git a/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py b/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py index 0bb3263496..ee8917fb0e 100644 --- a/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py @@ -15,8 +15,8 @@ from typing import List, Tuple import numpy as np import matplotlib.pyplot as plt -from qiskit.result import LocalReadoutMitigator from qiskit.result import marginal_distribution +from qiskit_experiments.data_processing import LocalReadoutMitigator from qiskit_experiments.framework import ExperimentData from qiskit_experiments.framework.matplotlib import get_non_gui_ax from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options diff --git a/qiskit_experiments/library/tomography/basis/pauli_basis.py b/qiskit_experiments/library/tomography/basis/pauli_basis.py index 1d908ef322..082da12333 100644 --- a/qiskit_experiments/library/tomography/basis/pauli_basis.py +++ b/qiskit_experiments/library/tomography/basis/pauli_basis.py @@ -16,9 +16,9 @@ import numpy as np from qiskit.circuit import QuantumCircuit from qiskit.circuit.library import HGate, XGate, ZGate, SGate, SdgGate -from qiskit.result import LocalReadoutMitigator from qiskit.quantum_info import DensityMatrix from qiskit.exceptions import QiskitError +from qiskit_experiments.data_processing import LocalReadoutMitigator from .local_basis import LocalMeasurementBasis, LocalPreparationBasis diff --git a/releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml b/releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml new file mode 100644 index 0000000000..61a06234a4 --- /dev/null +++ b/releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml @@ -0,0 +1,18 @@ +--- +features: + - | + New :class:`~.LocalReadoutMitigator` and + :class:`~.CorrelatedReadoutMitigator` classes have been added. These + classes were moved directly from Qiskit which deprecated them in Qiskit + 1.3. They provide utility methods for applying readout error mitigation and + integrate with the readout error mitigation experiments + :class:`~.LocalReadoutError` and :class:`~.CorrelatedReadoutError`. +upgrade: + - | + The readout error mitigation experiments :class:`~.LocalReadoutError` and + :class:`~.CorrelatedReadoutError` have been updated to generate instances + of the new :class:`~.LocalReadoutMitigator` and + :class:`~.CorrelatedReadoutMitigator` classes. The experiments should + continue to work as before, but any code that was using, for example, + ``isinstance()`` to check object type would need to be updated to check + against the Qiskit Experiments classes instead of the old Qiskit classes. diff --git a/test/data_processing/test_mitigators.py b/test/data_processing/test_mitigators.py new file mode 100644 index 0000000000..b75f7fde12 --- /dev/null +++ b/test/data_processing/test_mitigators.py @@ -0,0 +1,477 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +# pylint: disable=invalid-name + +"""Tests for error mitigation routines.""" + +import unittest +from collections import Counter +from test.base import QiskitExperimentsTestCase + +import numpy as np + +from qiskit import QiskitError, QuantumCircuit +from qiskit.quantum_info import Statevector +from qiskit.quantum_info.operators.predicates import matrix_equal +from qiskit.result import Counts +from qiskit.result.utils import marginal_counts +from qiskit.providers.fake_provider import Fake5QV1 + +from qiskit_experiments.data_processing import ( + CorrelatedReadoutMitigator, + LocalReadoutMitigator, +) +from qiskit_experiments.data_processing.mitigation import ( + counts_probability_vector, + expval_with_stddev, + stddev, + str2diag, +) + + +class TestReadoutMitigation(QiskitExperimentsTestCase): + """Tests for correlated and local readout mitigation.""" + + @classmethod + def setUpClass(cls): + super().setUpClass() + cls.rng = np.random.default_rng(42) + + @staticmethod + def compare_results(res1, res2): + """Compare the results between two runs""" + res1_total_shots = sum(res1.values()) + res2_total_shots = sum(res2.values()) + keys = set(res1.keys()).union(set(res2.keys())) + total = 0 + for key in keys: + val1 = res1.get(key, 0) / res1_total_shots + val2 = res2.get(key, 0) / res2_total_shots + total += abs(val1 - val2) ** 2 + return total + + @staticmethod + def mitigators(assignment_matrices, qubits=None): + """Generates the mitigators to test for given assignment matrices""" + full_assignment_matrix = assignment_matrices[0] + for m in assignment_matrices[1:]: + full_assignment_matrix = np.kron(full_assignment_matrix, m) + CRM = CorrelatedReadoutMitigator(full_assignment_matrix, qubits) + LRM = LocalReadoutMitigator(assignment_matrices, qubits) + mitigators = [CRM, LRM] + return mitigators + + @staticmethod + def simulate_circuit(circuit, assignment_matrix, num_qubits, shots=1024): + """Simulates the given circuit under the given readout noise""" + probs = Statevector.from_instruction(circuit).probabilities() + noisy_probs = assignment_matrix @ probs + labels = [bin(a)[2:].zfill(num_qubits) for a in range(2**num_qubits)] + results = TestReadoutMitigation.rng.choice(labels, size=shots, p=noisy_probs) + return Counts(dict(Counter(results))) + + @staticmethod + def ghz_3_circuit(): + """A 3-qubit circuit generating |000>+|111>""" + c = QuantumCircuit(3) + c.h(0) + c.cx(0, 1) + c.cx(1, 2) + return (c, "ghz_3_ciruit", 3) + + @staticmethod + def first_qubit_h_3_circuit(): + """A 3-qubit circuit generating |000>+|001>""" + c = QuantumCircuit(3) + c.h(0) + return (c, "first_qubit_h_3_circuit", 3) + + @staticmethod + def assignment_matrices(): + """A 3-qubit readout noise assignment matrices""" + return LocalReadoutMitigator(backend=Fake5QV1())._assignment_mats[0:3] + + @staticmethod + def counts_data(circuit, assignment_matrices, shots=1024): + """Generates count data for the noisy and noiseless versions of the circuit simulation""" + full_assignment_matrix = assignment_matrices[0] + for m in assignment_matrices[1:]: + full_assignment_matrix = np.kron(full_assignment_matrix, m) + num_qubits = len(assignment_matrices) + ideal_assignment_matrix = np.eye(2**num_qubits) + counts_ideal = TestReadoutMitigation.simulate_circuit( + circuit, ideal_assignment_matrix, num_qubits, shots + ) + counts_noise = TestReadoutMitigation.simulate_circuit( + circuit, full_assignment_matrix, num_qubits, shots + ) + probs_noise = {key: value / shots for key, value in counts_noise.items()} + return counts_ideal, counts_noise, probs_noise + + def test_mitigation_improvement(self): + """Test whether readout mitigation led to more accurate results""" + shots = 1024 + with self.assertWarns(DeprecationWarning): + # TODO self.assignment_matrices calls LocalReadoutMitigator, + # which only supports BackendV1 at the moment: + # https://github.com/Qiskit/qiskit/issues/12832 + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices) + circuit, circuit_name, num_qubits = self.ghz_3_circuit() + counts_ideal, counts_noise, probs_noise = self.counts_data( + circuit, assignment_matrices, shots + ) + unmitigated_error = self.compare_results(counts_ideal, counts_noise) + unmitigated_stddev = stddev(probs_noise, shots) + + for mitigator in mitigators: + mitigated_quasi_probs = mitigator.quasi_probabilities(counts_noise) + mitigated_probs = ( + mitigated_quasi_probs.nearest_probability_distribution().binary_probabilities( + num_bits=num_qubits + ) + ) + mitigated_error = self.compare_results(counts_ideal, mitigated_probs) + self.assertLess( + mitigated_error, + unmitigated_error * 0.8, + f"Mitigator {mitigator} did not improve circuit {circuit_name} measurements", + ) + mitigated_stddev_upper_bound = mitigated_quasi_probs._stddev_upper_bound + max_unmitigated_stddev = max(unmitigated_stddev.values()) + self.assertGreaterEqual( + mitigated_stddev_upper_bound, + max_unmitigated_stddev, + f"Mitigator {mitigator} on circuit {circuit_name} gave stddev upper bound " + f"{mitigated_stddev_upper_bound} while unmitigated stddev maximum is " + f"{max_unmitigated_stddev}", + ) + + def test_expectation_improvement(self): + """Test whether readout mitigation led to more accurate results + and that its standard deviation is increased""" + shots = 1024 + with self.assertWarns(DeprecationWarning): + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices) + num_qubits = len(assignment_matrices) + diagonals = [] + diagonals.append("IZ0") + diagonals.append("101") + diagonals.append("IZZ") + qubit_index = {i: i for i in range(num_qubits)} + circuit, circuit_name, num_qubits = self.ghz_3_circuit() + counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots) + probs_ideal, _ = counts_probability_vector(counts_ideal, qubit_index=qubit_index) + probs_noise, _ = counts_probability_vector(counts_noise, qubit_index=qubit_index) + for diagonal in diagonals: + if isinstance(diagonal, str): + diagonal = str2diag(diagonal) + unmitigated_expectation, unmitigated_stddev = expval_with_stddev( + diagonal, probs_noise, shots=counts_noise.shots() + ) + ideal_expectation = np.dot(probs_ideal, diagonal) + unmitigated_error = np.abs(ideal_expectation - unmitigated_expectation) + for mitigator in mitigators: + mitigated_expectation, mitigated_stddev = mitigator.expectation_value( + counts_noise, diagonal + ) + mitigated_error = np.abs(ideal_expectation - mitigated_expectation) + self.assertLess( + mitigated_error, + unmitigated_error, + f"Mitigator {mitigator} did not improve circuit {circuit_name} expectation " + f"computation for diagonal {diagonal} ideal: {ideal_expectation}, unmitigated:" + f" {unmitigated_expectation} mitigated: {mitigated_expectation}", + ) + self.assertGreaterEqual( + mitigated_stddev, + unmitigated_stddev, + f"Mitigator {mitigator} did not increase circuit {circuit_name} the" + f" standard deviation", + ) + + def test_clbits_parameter(self): + """Test whether the clbits parameter is handled correctly""" + shots = 10000 + with self.assertWarns(DeprecationWarning): + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices) + circuit, _, _ = self.first_qubit_h_3_circuit() + counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots) + counts_ideal_12 = marginal_counts(counts_ideal, [1, 2]) + counts_ideal_02 = marginal_counts(counts_ideal, [0, 2]) + + for mitigator in mitigators: + mitigated_probs_12 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[1, 2], clbits=[1, 2]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=2) + ) + mitigated_error = self.compare_results(counts_ideal_12, mitigated_probs_12) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly marginalize for qubits 1,2", + ) + + mitigated_probs_02 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[0, 2], clbits=[0, 2]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=2) + ) + mitigated_error = self.compare_results(counts_ideal_02, mitigated_probs_02) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly marginalize for qubits 0,2", + ) + + def test_qubits_parameter(self): + """Test whether the qubits parameter is handled correctly""" + shots = 10000 + with self.assertWarns(DeprecationWarning): + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices) + circuit, _, _ = self.first_qubit_h_3_circuit() + counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots) + counts_ideal_012 = counts_ideal + counts_ideal_210 = Counts({"000": counts_ideal["000"], "100": counts_ideal["001"]}) + counts_ideal_102 = Counts({"000": counts_ideal["000"], "010": counts_ideal["001"]}) + + for mitigator in mitigators: + mitigated_probs_012 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[0, 1, 2]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=3) + ) + mitigated_error = self.compare_results(counts_ideal_012, mitigated_probs_012) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly handle qubit order 0, 1, 2", + ) + + mitigated_probs_210 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[2, 1, 0]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=3) + ) + mitigated_error = self.compare_results(counts_ideal_210, mitigated_probs_210) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly handle qubit order 2, 1, 0", + ) + + mitigated_probs_102 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[1, 0, 2]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=3) + ) + mitigated_error = self.compare_results(counts_ideal_102, mitigated_probs_102) + self.assertLess( + mitigated_error, + 0.001, + "Mitigator {mitigator} did not correctly handle qubit order 1, 0, 2", + ) + + def test_repeated_qubits_parameter(self): + """Tests the order of mitigated qubits.""" + shots = 10000 + with self.assertWarns(DeprecationWarning): + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices, qubits=[0, 1, 2]) + circuit, _, _ = self.first_qubit_h_3_circuit() + counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots) + counts_ideal_012 = counts_ideal + counts_ideal_210 = Counts({"000": counts_ideal["000"], "100": counts_ideal["001"]}) + + for mitigator in mitigators: + mitigated_probs_210 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[2, 1, 0]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=3) + ) + mitigated_error = self.compare_results(counts_ideal_210, mitigated_probs_210) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly handle qubit order 2,1,0", + ) + + # checking qubit order 2,1,0 should not "overwrite" the default 0,1,2 + mitigated_probs_012 = ( + mitigator.quasi_probabilities(counts_noise) + .nearest_probability_distribution() + .binary_probabilities(num_bits=3) + ) + mitigated_error = self.compare_results(counts_ideal_012, mitigated_probs_012) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly handle qubit order 0,1,2 " + f"(the expected default)", + ) + + def test_qubits_subset_parameter(self): + """Tests mitigation on a subset of the initial set of qubits.""" + + shots = 10000 + with self.assertWarns(DeprecationWarning): + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices, qubits=[2, 4, 6]) + circuit, _, _ = self.first_qubit_h_3_circuit() + counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots) + counts_ideal_2 = marginal_counts(counts_ideal, [0]) + counts_ideal_6 = marginal_counts(counts_ideal, [2]) + + for mitigator in mitigators: + mitigated_probs_2 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[2]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=1) + ) + mitigated_error = self.compare_results(counts_ideal_2, mitigated_probs_2) + self.assertLess( + mitigated_error, + 0.001, + "Mitigator {mitigator} did not correctly handle qubit subset", + ) + + mitigated_probs_6 = ( + mitigator.quasi_probabilities(counts_noise, qubits=[6]) + .nearest_probability_distribution() + .binary_probabilities(num_bits=1) + ) + mitigated_error = self.compare_results(counts_ideal_6, mitigated_probs_6) + self.assertLess( + mitigated_error, + 0.001, + f"Mitigator {mitigator} did not correctly handle qubit subset", + ) + diagonal = str2diag("ZZ") + ideal_expectation = 0 + mitigated_expectation, _ = mitigator.expectation_value( + counts_noise, diagonal, qubits=[2, 6] + ) + mitigated_error = np.abs(ideal_expectation - mitigated_expectation) + self.assertLess( + mitigated_error, + 0.1, + f"Mitigator {mitigator} did not improve circuit expectation", + ) + + def test_from_backend(self): + """Test whether a local mitigator can be created directly from backend properties""" + with self.assertWarns(DeprecationWarning): + backend = Fake5QV1() + num_qubits = len(backend.properties().qubits) + probs = TestReadoutMitigation.rng.random((num_qubits, 2)) + for qubit_idx, qubit_prop in enumerate(backend.properties().qubits): + for prop in qubit_prop: + if prop.name == "prob_meas1_prep0": + prop.value = probs[qubit_idx][0] + if prop.name == "prob_meas0_prep1": + prop.value = probs[qubit_idx][1] + LRM_from_backend = LocalReadoutMitigator(backend=backend) + + mats = [] + for qubit_idx in range(num_qubits): + mat = np.array( + [ + [1 - probs[qubit_idx][0], probs[qubit_idx][1]], + [probs[qubit_idx][0], 1 - probs[qubit_idx][1]], + ] + ) + mats.append(mat) + LRM_from_matrices = LocalReadoutMitigator(assignment_matrices=mats) + self.assertTrue( + matrix_equal( + LRM_from_backend.assignment_matrix(), LRM_from_matrices.assignment_matrix() + ) + ) + + def test_error_handling(self): + """Test that the assignment matrices are valid.""" + bad_matrix_A = np.array([[-0.3, 1], [1.3, 0]]) # negative indices + bad_matrix_B = np.array([[0.2, 1], [0.7, 0]]) # columns not summing to 1 + good_matrix_A = np.array([[0.2, 1], [0.8, 0]]) + for bad_matrix in [bad_matrix_A, bad_matrix_B]: + with self.assertRaises(QiskitError) as cm: + CorrelatedReadoutMitigator(bad_matrix) + self.assertEqual( + cm.exception.message, + "Assignment matrix columns must be valid probability distributions", + ) + + with self.assertRaises(QiskitError) as cm: + amats = [good_matrix_A, bad_matrix_A] + LocalReadoutMitigator(amats) + self.assertEqual( + cm.exception.message, + "Assignment matrix columns must be valid probability distributions", + ) + + with self.assertRaises(QiskitError) as cm: + amats = [bad_matrix_B, good_matrix_A] + LocalReadoutMitigator(amats) + self.assertEqual( + cm.exception.message, + "Assignment matrix columns must be valid probability distributions", + ) + + def test_expectation_value_endian(self): + """Test that endian for expval is little.""" + with self.assertWarns(DeprecationWarning): + assignment_matrices = self.assignment_matrices() + mitigators = self.mitigators(assignment_matrices) + counts = Counts({"10": 3, "11": 24, "00": 74, "01": 923}) + for mitigator in mitigators: + expval, _ = mitigator.expectation_value(counts, diagonal="IZ", qubits=[0, 1]) + self.assertAlmostEqual(expval, -1.0, places=0) + + def test_quasi_probabilities_shots_passing(self): + """Test output of LocalReadoutMitigator.quasi_probabilities + + We require the number of shots to be set in the output. + """ + mitigator = LocalReadoutMitigator([np.array([[0.9, 0.1], [0.1, 0.9]])], qubits=[0]) + counts = Counts({"10": 3, "11": 24, "00": 74, "01": 923}) + quasi_dist = mitigator.quasi_probabilities(counts) + self.assertEqual(quasi_dist.shots, sum(counts.values())) + + # custom number of shots + quasi_dist = mitigator.quasi_probabilities(counts, shots=1025) + self.assertEqual(quasi_dist.shots, 1025) + + +class TestLocalReadoutMitigation(QiskitExperimentsTestCase): + """Tests specific to the local readout mitigator""" + + def test_assignment_matrix(self): + """Tests that the local mitigator generates the full assignment matrix correctly""" + qubits = [7, 2, 3] + with self.assertWarns(DeprecationWarning): + backend = Fake5QV1() + assignment_matrices = LocalReadoutMitigator(backend=backend)._assignment_mats[0:3] + expected_assignment_matrix = np.kron( + np.kron(assignment_matrices[2], assignment_matrices[1]), assignment_matrices[0] + ) + expected_mitigation_matrix = np.linalg.inv(expected_assignment_matrix) + LRM = LocalReadoutMitigator(assignment_matrices, qubits) + self.assertTrue(matrix_equal(expected_mitigation_matrix, LRM.mitigation_matrix())) + self.assertTrue(matrix_equal(expected_assignment_matrix, LRM.assignment_matrix())) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/library/tomography/test_process_tomography.py b/test/library/tomography/test_process_tomography.py index 18ae55dac8..bdc0fa3c0a 100644 --- a/test/library/tomography/test_process_tomography.py +++ b/test/library/tomography/test_process_tomography.py @@ -22,11 +22,11 @@ import qiskit.quantum_info as qi from qiskit import QuantumCircuit from qiskit.circuit.library import XGate, CXGate -from qiskit.result import LocalReadoutMitigator from qiskit_aer import AerSimulator from qiskit_aer.noise import NoiseModel +from qiskit_experiments.data_processing import LocalReadoutMitigator from qiskit_experiments.database_service import ExperimentEntryNotFound from qiskit_experiments.library import ProcessTomography, MitigatedProcessTomography from qiskit_experiments.library.tomography import ProcessTomographyAnalysis, basis diff --git a/test/library/tomography/test_state_tomography.py b/test/library/tomography/test_state_tomography.py index b0ec4b52c8..d4c3cb2ec4 100644 --- a/test/library/tomography/test_state_tomography.py +++ b/test/library/tomography/test_state_tomography.py @@ -23,11 +23,11 @@ import qiskit.quantum_info as qi from qiskit import QuantumCircuit from qiskit.circuit.library import XGate -from qiskit.result import LocalReadoutMitigator from qiskit_aer import AerSimulator from qiskit_aer.noise import NoiseModel +from qiskit_experiments.data_processing import LocalReadoutMitigator from qiskit_experiments.database_service import ExperimentEntryNotFound from qiskit_experiments.library import StateTomography, MitigatedStateTomography from qiskit_experiments.library.tomography import StateTomographyAnalysis, basis From 519bb534c134ca9a251f4f92f6abd39f56017341 Mon Sep 17 00:00:00 2001 From: Luciano Bello Date: Fri, 15 Nov 2024 17:57:40 +0100 Subject: [PATCH 2/2] Move Provider import only under `TYPE_CHECKING` conditions (#1489) ### Summary Qiskit 1.1 deprecated the Provider abstract class to be removed in 2.0 (planned for 1Q2025). In those cases where the class is used for only typechecking purposes, the import should only happen when `TYPE_CHECKING`. For the rest of the cases, I opened https://github.com/qiskit-community/qiskit-experiments/issues/1488 --- qiskit_experiments/framework/experiment_data.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/qiskit_experiments/framework/experiment_data.py b/qiskit_experiments/framework/experiment_data.py index db71273ffd..a4cd953967 100644 --- a/qiskit_experiments/framework/experiment_data.py +++ b/qiskit_experiments/framework/experiment_data.py @@ -36,7 +36,7 @@ from qiskit.result import Result from qiskit.providers.jobstatus import JobStatus, JOB_FINAL_STATES from qiskit.exceptions import QiskitError -from qiskit.providers import Job, Backend, Provider +from qiskit.providers import Job, Backend from qiskit.utils.deprecation import deprecate_arg from qiskit.primitives import BitArray, SamplerPubResult, BasePrimitiveJob @@ -76,6 +76,7 @@ # `TYPE_CHECKING` means that the import will never be resolved by an actual # interpreter, only static analysis. from . import BaseExperiment + from qiskit.providers import Provider LOG = logging.getLogger(__name__)