From 2f0d192fa7848b804dfc02acc101a260083dbe85 Mon Sep 17 00:00:00 2001 From: Will Shanks Date: Tue, 5 Nov 2024 17:27:13 -0500 Subject: [PATCH] Adopt qiskit.result.mitigation into qiskit_experiments.data_processing --- .../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 | 159 ++++++ .../correlated_readout_error_analysis.py | 2 +- .../analysis/local_readout_error_analysis.py | 2 +- .../library/tomography/basis/pauli_basis.py | 2 +- test/data_processing/test_mitigators.py | 477 ++++++++++++++++++ .../tomography/test_process_tomography.py | 2 +- .../tomography/test_state_tomography.py | 2 +- 12 files changed, 1346 insertions(+), 5 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 test/data_processing/test_mitigators.py 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..ebb97d078b --- /dev/null +++ b/qiskit_experiments/data_processing/mitigation/utils.py @@ -0,0 +1,159 @@ +# 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 +""" + +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 marganalize_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 = marganalize_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/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