diff --git a/README.md b/README.md index e963d03e0..f359758f6 100644 --- a/README.md +++ b/README.md @@ -26,9 +26,13 @@ visualization functions for fitters you'll need to install matplotlib. You can do this with `pip install matplotlib` or when you install ignis with `pip install qiskit-ignis[visualization]`. If you're going to use a cvx fitter for running tomogography you'll need to install cvxpy. You can do this with -`pip install cvxpy` or when you install install ignis with -`pip install qiskit-ignis[cvx]`. If you want to install both when you install -ignis you can run `pip install qiskit-ignis[visualization,cvx]`. +`pip install cvxpy` or when you install ignis with +`pip install qiskit-ignis[cvx]`. When performing expectation value measurement +error mitigation using the CTMP method performance can be improved using +just-in-time compiling if Numbda is installed. You can do this with +`pip install numba` or when you install ignis with +`pip install qiskit-ignis[jit]`. If you want to install all extra requirements +when you install ignis you can run `pip install qiskit-ignis[visualization,cvx,jit]`. ## Creating your first quantum experiment with Qiskit Ignis Now that you have Qiskit Ignis installed, you can start creating experiments, to reveal information about the device quality. Here is a basic example: diff --git a/qiskit/ignis/mitigation/__init__.py b/qiskit/ignis/mitigation/__init__.py index ec8f2348f..dcf756ba2 100644 --- a/qiskit/ignis/mitigation/__init__.py +++ b/qiskit/ignis/mitigation/__init__.py @@ -36,7 +36,30 @@ CompleteMeasFitter TensoredMeasFitter +Expectation Value Measurement +============================= + +The following classes allow mitigation of measurement errors when computing +expectation values of diagonal operators from counts. + +.. autosummary:: + :toctree: ../stubs/ + + expectation_value + expval_meas_mitigator_circuits + ExpvalMeasMitigatorFitter + CTMPExpvalMeasMitigator + CompleteExpvalMeasMitigator + TensoredExpvalMeasMitigator """ + from .measurement import (complete_meas_cal, tensored_meas_cal, MeasurementFilter, TensoredFilter, CompleteMeasFitter, TensoredMeasFitter) + +from .expval import (expectation_value, + expval_meas_mitigator_circuits, + ExpvalMeasMitigatorFitter, + CompleteExpvalMeasMitigator, + TensoredExpvalMeasMitigator, + CTMPExpvalMeasMitigator) diff --git a/qiskit/ignis/mitigation/expval/__init__.py b/qiskit/ignis/mitigation/expval/__init__.py new file mode 100644 index 000000000..b735ff078 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/__init__.py @@ -0,0 +1,23 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. +""" +Expectation value measurement error mitigation module +""" + +from .utils import expectation_value +from .circuits import expval_meas_mitigator_circuits +from .fitter import ExpvalMeasMitigatorFitter +from .complete_mitigator import CompleteExpvalMeasMitigator +from .tensored_mitigator import TensoredExpvalMeasMitigator +from .ctmp_mitigator import CTMPExpvalMeasMitigator diff --git a/qiskit/ignis/mitigation/expval/base_meas_mitigator.py b/qiskit/ignis/mitigation/expval/base_meas_mitigator.py new file mode 100644 index 000000000..c6b445b3f --- /dev/null +++ b/qiskit/ignis/mitigation/expval/base_meas_mitigator.py @@ -0,0 +1,271 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. +""" +Full A-matrix measurement migitation generator. +""" + +from abc import ABC, abstractmethod +from typing import Optional, List, Dict, Tuple, Union +import numpy as np + +try: + import matplotlib.pyplot as plt + _HAS_MATPLOTLIB = True +except ImportError: + _HAS_MATPLOTLIB = False + + +class BaseExpvalMeasMitigator(ABC): + """Base measurement error mitigator class.""" + + @abstractmethod + def expectation_value(self, + counts: Dict, + diagonal: Optional[ + Union[np.ndarray, List[complex], List[float]]] = None, + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, + ) -> Tuple[float, float]: + r"""Compute a measurement error mitigated expectation value. + + 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: + counts: counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the 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. + + Returns: + (float, float): the expectation value and 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)``. + """ + + @abstractmethod + def mitigation_matrix(self, qubits: List[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. + + Returns: + np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`. + """ + + @abstractmethod + 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. + """ + + @abstractmethod + def _compute_gamma(self, qubits: Optional[List[int]] = None) -> float: + """Compute gamma for N-qubit mitigation.""" + + def assignment_fidelity(self, qubits: Optional[List[int]] = None) -> float: + r"""Return the measurement assignment fidelity on the specified qubits. + + The assignment fidelity on N-qubits is defined as + :math:`\sum_{x\in\{0, 1\}^n} P(x|x) / 2^n`, where + :math:`P(x|x) = \rangle x|A|x\langle`, and :math:`A` is the + :meth:`assignment_matrix`. + + Args: + qubits: Optional, qubits being measured for operator expval. + + Returns: + float: the assignment fidelity. + """ + return self.assignment_matrix(qubits=qubits).diagonal().mean() + + def required_shots(self, delta: float, qubits: Optional[List[int]] = None) -> int: + r"""Return the number of shots required for expectation value estimation. + + This is the number of shots required so that + :math:`|\langle O \rangle_{est} - \langle O \rangle_{true}| < \delta` + with high probability (at least 2/3) and is given by + :math:`4\delta^2 \Gamma^2` where :math:`\Gamma^2` is the + :meth:`mitigation_overhead`. + + Args: + delta: Error tolerance for expectation value estimator. + qubits: Optional, qubits being measured for operator expval. + + Returns: + int: the required shots. + """ + gamma = self._compute_gamma(qubits=qubits) + return int(np.ceil((4 * gamma ** 2) / (delta ** 2))) + + def mitigation_overhead(self, qubits: Optional[List[int]] = None) -> int: + """Return the mitigation overhead for expectation value estimation. + + This is the multiplicative factor of extra shots required for + estimating a mitigated expectation value with the same accuracy as + an unmitigated expectation value. + + Args: + qubits: Optional, qubits being measured for operator expval. + + Returns: + int: the mitigation overhead factor. + """ + gamma = self._compute_gamma(qubits=qubits) + return int(np.ceil(gamma ** 2)) + + def stddev_upper_bound(self, shots: int = 1, qubits: Optional[List[int]] = None) -> float: + """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 / np.sqrt(shots) + + def plot_assignment_matrix(self, + qubits: Optional[List[int]] = None, + ax: Optional[plt.axes] = None) -> plt.axes: + """Matrix plot of the readout error assignment matrix. + + Args: + qubits: Optional, qubits being measured for operator expval. + ax: Optional. Axes object to add plot to. + + Returns: + plt.axes: the figure axes object. + + Raises: + ImportError: if matplotlib is not installed. + """ + if not _HAS_MATPLOTLIB: + raise ImportError( + 'Plotting functions require the optional matplotlib package.' + ' Install with `pip install matplotlib`.') + + mat = self.assignment_matrix(qubits) + if ax is None: + figure = plt.figure() + ax = figure.gca() + axim = ax.matshow(mat, cmap=plt.cm.binary, clim=[0, 1]) + ax.figure.colorbar(axim) + ax = self._plot_axis(mat, ax=ax) + return ax + + def plot_mitigation_matrix(self, + qubits: Optional[List[int]] = None, + ax: Optional[plt.axes] = None) -> plt.axes: + """Matrix plot of the readout error mitigation matrix. + + Args: + qubits: Optional, qubits being measured for operator expval. + ax: Optional. Axes object to add plot to. + + Returns: + plt.axes: the figure axes object. + + Raises: + ImportError: if matplotlib is not installed. + """ + if not _HAS_MATPLOTLIB: + raise ImportError( + 'Plotting functions require the optional matplotlib package.' + ' Install with `pip install matplotlib`.') + + # Matrix parameters + mat = self.mitigation_matrix(qubits) + mat_max = np.max(mat) + mat_min = np.min(mat) + lim = max(abs(mat_max), abs(mat_min), 1) + + if ax is None: + figure = plt.figure() + ax = figure.gca() + axim = ax.matshow(mat, cmap=plt.cm.RdGy, clim=[-lim, lim]) + ax.figure.colorbar(axim, values=np.linspace(mat_min, mat_max, 100)) + ax = self._plot_axis(mat, ax=ax) + return ax + + @staticmethod + 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) + + @staticmethod + def _int_to_bitstring(i, num_qubits=None): + """Convert an integer to a bitstring.""" + label = bin(i)[2:] + if num_qubits: + pad = num_qubits - len(label) + if pad > 0: + label = pad * '0' + label + return label + + @staticmethod + def _plot_axis(mat: np.ndarray, ax: plt.axes) -> plt.axes: + """Helper function for setting up axes for plots. + + Args: + mat: the N-qubit matrix to plot. + ax: Optional. Axes object to add plot to. + + Returns: + plt.axes: the figure object and axes object. + """ + dim = len(mat) + num_qubits = int(np.log2(dim)) + bit_labels = [BaseExpvalMeasMitigator._int_to_bitstring(i, num_qubits) for i in range(dim)] + + ax.set_xticks(np.arange(dim)) + ax.set_yticks(np.arange(dim)) + ax.set_xticklabels(bit_labels, rotation=90) + ax.set_yticklabels(bit_labels) + ax.set_xlabel('Prepared State') + ax.xaxis.set_label_position('top') + ax.set_ylabel('Measured State') + return ax diff --git a/qiskit/ignis/mitigation/expval/circuits.py b/qiskit/ignis/mitigation/expval/circuits.py new file mode 100644 index 000000000..777382108 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/circuits.py @@ -0,0 +1,196 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. +""" +Expectation value measurement error migitation generator. +""" + +from typing import Optional, Tuple, List, Dict + +from qiskit import QuantumCircuit +from qiskit.exceptions import QiskitError + + +def expval_meas_mitigator_circuits(num_qubits: int, + method: Optional[str] = 'CTMP', + labels: Optional[List[str]] = None) -> Tuple[ + List[QuantumCircuit], List[Dict[str, any]] + ]: + """Generate measurement error mitigator circuits and metadata. + + Use the :class:`~qiskit.ignis.mitigation.ExpvalMeasMitigatorFitter` + class to fit the execution results to construct a calibrated expectation + value measurement error mitigator. + + Args: + num_qubits: the number of qubits to calibrate. + method: the mitigation method ``'complete'``, ``'tensored'``, or ``'CTMP'``. + labels: Optional, custom labels to run for calibration. If None + the method will determine the default label values. + + Returns: + tuple: (circuits, metadata) the measurement error characterization + circuits, and metadata for the fitter. + + Mitigation Method: + * The ``'complete'`` method will generate all :math:`2^n` computational + basis states measurement circuits and fitting will return a + :class:`~qiskit.ignis.mitigation.CompleteExpvalMeasMitigator`. This + method should only be used for small numbers of qubits. + * The ``'tensored'`` method will generate two input state circuits of + the all 0 and all 1 states on number of qubits unless custom labels + are specified. Ftting will return a + :class:`~qiskit.ignis.mitigation.TensoredExpvalMeasMitigator`. This + method assumes measurement errors are uncorrelated between qubits. + * The ``'CTMP'`` method will generate :math:`n+2` input state circuits + unless custom labels are specified. The default input states are + the all 0 state, the all 1 state, and the :math:`n` state with a + single qubit in the 1 state and all others in the 0 state. + Ftting will return a + :class:`~qiskit.ignis.mitigation.CTMPExpvalMeasMitigator`. + + Example: + + The following example shows calibrating a 5-qubit expectation value + measurement error mitigator using the ``'tensored'`` method. + + .. jupyter-execute:: + + from qiskit import execute + from qiskit.test.mock import FakeVigo + import qiskit.ignis.mitigation as mit + + backend = FakeVigo() + num_qubits = backend.configuration().num_qubits + + # Generate calibration circuits + circuits, metadata = mit.expval_meas_mitigator_circuits( + num_qubits, method='tensored') + result = execute(circuits, backend, shots=8192).result() + + # Fit mitigator + mitigator = mit.ExpvalMeasMitigatorFitter(result, metadata).fit() + + # Plot fitted N-qubit assignment matrix + mitigator.plot_assignment_matrix() + + The following shows how to use the above mitigator to apply measurement + error mitigation to expectation value computations + + .. jupyter-execute:: + + from qiskit import QuantumCircuit + + # Test Circuit with expectation value -1. + qc = QuantumCircuit(num_qubits) + qc.x(range(num_qubits)) + qc.measure_all() + + # Execute + shots = 8192 + seed_simulator = 1999 + result = execute(qc, backend, shots=8192, seed_simulator=1999).result() + counts = result.get_counts(0) + + # Expectation value of Z^N without mitigation + expval_nomit, error_nomit = mit.expectation_value(counts) + print('Expval (no mitigation): {:.2f} \u00B1 {:.2f}'.format( + expval_nomit, error_nomit)) + + # Expectation value of Z^N with mitigation + expval_mit, error_mit = mit.expectation_value(counts, + meas_mitigator=mitigator) + print('Expval (with mitigation): {:.2f} \u00B1 {:.2f}'.format( + expval_mit, error_mit)) + + """ + generator = ExpvalMeasMitigatorCircuits(num_qubits, method, labels) + return generator.generate_circuits() + + +class ExpvalMeasMitigatorCircuits: + """Expecation value measurement error mitigator calibration circuits.""" + + # pylint: disable=arguments-differ + def __init__(self, + num_qubits: int, + method: str = 'CTMP', + labels: Optional[List[str]] = None): + """Initialize measurement mitigator calibration generator. + + Args: + num_qubits: the number of qubits to calibrate. + method: the mitigation method 'complete', 'tensored', or 'CTMP'. + labels: custom labels to run for calibration. + """ + self._num_qubits = num_qubits + self._circuits = [] + self._metadata = [] + if labels is None: + labels = self._method_labels(method) + for label in labels: + self._metadata.append({ + 'experiment': 'meas_mit', + 'cal': label, + 'method': method, + }) + self._circuits.append(self._calibration_circuit(num_qubits, label)) + + def _method_labels(self, method): + """Generate labels for initializing via a standard method.""" + + if method == 'tensored': + return [self._num_qubits * '0', self._num_qubits * '1'] + + if method in ['CTMP', 'ctmp']: + labels = [self._num_qubits * '0', self._num_qubits * '1'] + for i in range(self._num_qubits): + labels.append(((self._num_qubits - i - 1) * '0') + '1' + + (i * '0')) + return labels + + if method == 'complete': + labels = [] + for i in range(2**self._num_qubits): + bits = bin(i)[2:] + label = (self._num_qubits - len(bits)) * '0' + bits + labels.append(label) + return labels + + raise QiskitError("Unrecognized method {}".format(method)) + + def generate_circuits(self) -> Tuple[List[QuantumCircuit], List[dict]]: + """Return experiment payload data. +​ + Circuits is a list of circuits for the experiments, metadata is a list of metadata + for the experiment that is required by the fitter to interpreting results. +​ + Returns: + tuple: circuits, metadata + """ + return self._circuits, self._metadata + + @staticmethod + def _calibration_circuit(num_qubits: int, label: str) -> QuantumCircuit: + """Return a calibration circuit. + + This is an N-qubit circuit where N is the length of the label. + The circuit consists of X-gates on qubits with label bits equal to 1, + and measurements of all qubits. + """ + circ = QuantumCircuit(num_qubits, name='meas_mit_cal_' + label) + for i, val in enumerate(reversed(label)): + if val == '1': + circ.x(i) + circ.measure_all() + return circ diff --git a/qiskit/ignis/mitigation/expval/complete_mitigator.py b/qiskit/ignis/mitigation/expval/complete_mitigator.py new file mode 100644 index 000000000..5d6b67ad9 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/complete_mitigator.py @@ -0,0 +1,180 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +Full-matrix measurement error mitigation generator. +""" +from typing import Optional, List, Dict, Tuple +import numpy as np + +from qiskit.exceptions import QiskitError +from .utils import (counts_probability_vector, _expval_with_stddev) +from .base_meas_mitigator import BaseExpvalMeasMitigator + + +class CompleteExpvalMeasMitigator(BaseExpvalMeasMitigator): + """N-qubit measurement error mitigator. + + This class can be used with the + :func:`qiskit.ignis.mitigation.expectation_value` function to apply + measurement error mitigation of general N-qubit measurement errors + when calculating expectation values from counts. Expectation values + can also be computed directly using the :meth:`expectation_value` method. + + For measurement mitigation to be applied the mitigator should be + calibrated using the + :func:`qiskit.ignis.mitigation.expval_meas_mitigator_circuits` function + and :class:`qiskit.ignis.mitigation.ExpvalMeasMitigatorFitter` class with + the ``'complete'`` mitigation method. + """ + + def __init__(self, amat: np.ndarray): + """Initialize a TensorMeasurementMitigator + + Args: + amat (np.array): readout error assignment matrix. + """ + self._num_qubits = int(np.log2(amat.shape[0])) + self._assignment_mat = amat + self._mitigation_mats = {} + + def expectation_value(self, + counts: Dict, + diagonal: Optional[np.ndarray] = None, + qubits: Optional[List[int]] = None, + clbits: Optional[List[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: + counts: counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the 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. + + Returns: + (float, float): the expectation value and standard deviation. + + Raises: + QiskitError: if input arguments are invalid. + + 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)``. + """ + # Get probability vector + probs, shots = counts_probability_vector( + counts, clbits=clbits, qubits=qubits, return_shots=True) + num_qubits = int(np.log2(probs.shape[0])) + + # Get qubit mitigation matrix and mitigate probs + if qubits is None: + qubits = tuple(range(num_qubits)) + if len(qubits) != num_qubits: + raise QiskitError("Num qubits does not match number of clbits.") + mit_mat = self.mitigation_matrix(qubits) + + # Get operator coeffs + if diagonal is None: + diagonal = self._z_diagonal(2 ** num_qubits) + # Apply transpose of mitigation matrix + coeffs = mit_mat.T.dot(diagonal) + + return _expval_with_stddev(coeffs, probs, shots) + + def mitigation_matrix(self, qubits: List[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. + + Returns: + np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`. + """ + if qubits is None: + qubits = tuple(range(self._num_qubits)) + else: + 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 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: + return self._assignment_mat + + if isinstance(qubits, int): + qubits = [qubits] + + # Compute marginal matrix + axis = tuple([self._num_qubits - 1 - i for i in set( + range(self._num_qubits)).difference(qubits)]) + 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(qubits)]): + 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, qubits=None): + """Compute gamma for N-qubit mitigation""" + mitmat = self.mitigation_matrix(qubits=qubits) + return np.max(np.sum(np.abs(mitmat), axis=0)) diff --git a/qiskit/ignis/mitigation/expval/ctmp_fitter.py b/qiskit/ignis/mitigation/expval/ctmp_fitter.py new file mode 100644 index 000000000..50cce64ee --- /dev/null +++ b/qiskit/ignis/mitigation/expval/ctmp_fitter.py @@ -0,0 +1,162 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +Full-matrix measurement error mitigation generator. +""" +import logging +from typing import List, Dict +import numpy as np +import scipy.linalg as la + +from qiskit.exceptions import QiskitError +from .ctmp_mitigator import CTMPExpvalMeasMitigator +from .ctmp_generator_set import Generator, standard_generator_set +from .utils import assignment_matrix + +logger = logging.getLogger(__name__) + + +def fit_ctmp_meas_mitigator(cal_data: Dict[int, Dict[int, int]], + num_qubits: int, + generators: List[Generator] = None) -> CTMPExpvalMeasMitigator: + """Return FullMeasureErrorMitigator from result data. + + Args: + cal_data: calibration dataset. + num_qubits: the number of qubits for the calibation dataset. + generators: Optional, input generator set. + + Returns: + Measurement error mitigator object. + + Raises: + QiskitError: if input arguments are invalid. + """ + if not isinstance(num_qubits, int): + raise QiskitError('Number of qubits must be an int') + if generators is None: + generators = standard_generator_set(num_qubits) + + gen_mat_dict = {} + for gen in generators + _supplementary_generators(generators): + if len(gen[2]) > 1: + mat = _local_g_matrix(gen, cal_data, num_qubits) + gen_mat_dict[gen] = mat + + # Compute rates for generators + rates = [_get_ctmp_error_rate(gen, gen_mat_dict, num_qubits) for gen in generators] + return CTMPExpvalMeasMitigator(generators, rates) + + +# Utility functions used for fitting (Should be moved to Fitter class) + +def _ctmp_err_rate_1_q(a: str, b: str, j: int, + g_mat_dict: Dict[Generator, np.array], + num_qubits: int) -> float: + """Compute the 1q error rate for a given generator.""" + # pylint: disable=invalid-name + rate_list = [] + if a == '0' and b == '1': + g1 = ('00', '10') + g2 = ('01', '11') + elif a == '1' and b == '0': + g1 = ('10', '00') + g2 = ('11', '01') + else: + raise ValueError('Invalid a,b encountered...') + for k in range(num_qubits): + if k == j: + continue + c = (j, k) + for g_strs in [g1, g2]: + gen = g_strs + (c,) + r = _ctmp_err_rate_2_q(gen, g_mat_dict) + rate_list.append(r) + if len(rate_list) != 2 * (num_qubits - 1): + raise ValueError('Rate list has wrong number of elements') + rate = np.mean(rate_list) + return rate + + +def _ctmp_err_rate_2_q(gen, g_mat_dict) -> float: + """Compute the 2 qubit error rate for a given generator.""" + # pylint: disable=invalid-name + g_mat = g_mat_dict[gen] + b, a, _ = gen + r = g_mat[int(b, 2), int(a, 2)] + return r + + +def _get_ctmp_error_rate(gen: Generator, + g_mat_dict: Dict[Generator, np.array], + num_qubits: int) -> float: + """Compute the error rate r_i for generator G_i. + + Args: + gen (Generator): Generator to calibrate. + g_mat_dict (Dict[Generator, np.array]): Dictionary of local G(j,k) matrices. + num_qubits (int): number of qubits. + + Returns: + float: The coefficient r_i for generator G_i. + + Raises: + ValueError: The provided generator is not already in the set of generators. + """ + # pylint: disable=invalid-name + b, a, c = gen + if len(b) == 1: + rate = _ctmp_err_rate_1_q(a, b, c[0], g_mat_dict, num_qubits) + elif len(b) == 2: + rate = _ctmp_err_rate_2_q(gen, g_mat_dict) + return rate + + +def _local_g_matrix(gen: Generator, + cal_data: Dict[int, Dict[int, int]], + num_qubits: int) -> np.array: + """Computes the G(j,k) matrix in the basis [00, 01, 10, 11].""" + # pylint: disable=invalid-name + _, _, c = gen + j, k = c + a_mat = assignment_matrix(cal_data, num_qubits, [j, k]) + g = la.logm(a_mat) + if np.linalg.norm(np.imag(g)) > 1e-3: + raise QiskitError('Encountered complex entries in G_i={}'.format(g)) + g = np.real(g) + for i in range(4): + for j in range(4): + if i != j: + if g[i, j] < 0: + g[i, j] = 0 + return g + + +def _supplementary_generators(gen_list: List[Generator]) -> List[Generator]: + """Supplementary generators needed to run 1q calibrations. + + Args: + gen_list (List[Generator]): List of generators. + + Returns: + List[Generator]: List of additional generators needed. + """ + pairs = {tuple(gen[2]) for gen in gen_list} + supp_gens = [] + for tup in pairs: + supp_gens.append(('10', '00', tup)) + supp_gens.append(('00', '10', tup)) + supp_gens.append(('11', '01', tup)) + supp_gens.append(('01', '11', tup)) + return supp_gens diff --git a/qiskit/ignis/mitigation/expval/ctmp_generator_set.py b/qiskit/ignis/mitigation/expval/ctmp_generator_set.py new file mode 100644 index 000000000..af9d8174b --- /dev/null +++ b/qiskit/ignis/mitigation/expval/ctmp_generator_set.py @@ -0,0 +1,91 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +"""Perform CTMP calibration for error mitigation. +""" + +from itertools import combinations +from typing import List, Tuple + +Generator = Tuple[str, str, Tuple[int]] + + +def standard_generator_set(num_qubits: int, pairs=None) -> List[Generator]: + """Construct this generator set given a number of qubits. + + Set of generators on 1 and 2 qubits. Corresponds to the following readout errors: + `0 -> 1` + `1 -> 0` + `01 -> 10` + `11 -> 00` + `00 -> 11` + + Args: + num_qubits (int): Number of qubits to calibrate. + pairs (Optional[List[Tuple[int, int]]]): The pairs of qubits to compute generators for. + If None, then use all pairs. Defaults to None. + + Returns: + StandardGeneratorSet: The generator set. + + Raises: + ValueError: Either the number of qubits is negative or not an int. + ValueError: The number of generators does not match what it should be by counting. + """ + if not isinstance(num_qubits, int): + raise ValueError('num_qubits needs to be an int') + if num_qubits <= 0: + raise ValueError('Need num_qubits at least 1') + generators = [] + generators += single_qubit_bitstrings(num_qubits) + if num_qubits > 1: + generators += two_qubit_bitstrings_symmetric(num_qubits, pairs=pairs) + generators += two_qubit_bitstrings_asymmetric(num_qubits, pairs=pairs) + if len(generators) != 2 * num_qubits ** 2: + raise ValueError('Incorrect length of generators. {} != 2n^2.'.format(len(generators))) + return generators + + +def single_qubit_bitstrings(num_qubits: int) -> List[Generator]: + """Returns a list of tuples `[(C_1, b_1, a_1), (C_2, b_2, a_2), ...]` that represent + the generators . + """ + res = [('1', '0', (i,)) for i in range(num_qubits)] + res += [('0', '1', (i,)) for i in range(num_qubits)] + if len(res) != 2 * num_qubits: + raise ValueError('Should have gotten 2n qubits, got {}'.format(len(res))) + return res + + +def two_qubit_bitstrings_symmetric(num_qubits: int, pairs=None) -> List[Generator]: + """Return the 11->00 and 00->11 generators on a given number of qubits. + """ + if pairs is None: + pairs = list(combinations(range(num_qubits), r=2)) + res = [('11', '00', (i, j)) for i, j in pairs if i < j] + res += [('00', '11', (i, j)) for i, j in pairs if i < j] + if len(res) != num_qubits * (num_qubits - 1): + raise ValueError('Should have gotten n(n-1) qubits, got {}'.format(len(res))) + return res + + +def two_qubit_bitstrings_asymmetric(num_qubits: int, pairs=None) -> List[Generator]: + """Return the 01->10 generators on a given number of qubits. + """ + if pairs is None: + pairs = list(combinations(range(num_qubits), r=2)) + res = [('10', '01', (i, j)) for i, j in pairs] + res += [('10', '01', (j, i)) for i, j in pairs] + if len(res) != num_qubits * (num_qubits - 1): + raise ValueError('Should have gotten n(n-1) qubits, got {}'.format(len(res))) + return res diff --git a/qiskit/ignis/mitigation/expval/ctmp_mitigator.py b/qiskit/ignis/mitigation/expval/ctmp_mitigator.py new file mode 100644 index 000000000..a50075b61 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/ctmp_mitigator.py @@ -0,0 +1,425 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +CTMP expectation value measurement error mitigator. +""" +import logging +from typing import Dict, Optional, List, Tuple + +import numpy as np +import scipy.linalg as la +import scipy.sparse as sps + +from qiskit.exceptions import QiskitError +from qiskit.ignis.numba import jit_fallback + +from .base_meas_mitigator import BaseExpvalMeasMitigator +from .utils import counts_probability_vector +from .ctmp_generator_set import Generator + +logger = logging.getLogger(__name__) + + +class CTMPExpvalMeasMitigator(BaseExpvalMeasMitigator): + """N-qubit CTMP measurement error mitigator. + + This class can be used with the + :func:`qiskit.ignis.mitigation.expectation_value` function to apply + measurement error mitigation of N-qubit measurement errors caused by + one and two-body error generators. Expectation values + can also be computed directly using the :meth:`expectation_value` method. + + For measurement mitigation to be applied the mitigator should be + calibrated using the + :func:`qiskit.ignis.mitigation.expval_meas_mitigator_circuits` function + and :class:`qiskit.ignis.mitigation.ExpvalMeasMitigatorFitter` class with + the ``'CTMP'`` mitigation method. + """ + def __init__(self, + generators: List[Generator], + rates: List[float], + num_qubits: Optional[int] = None, + seed: Optional = None): + """Initialize a TensorMeasurementMitigator""" + if num_qubits is None: + self._num_qubits = 1 + max([max([max(gen[2]) for gen in generators])]) + else: + self._num_qubits = num_qubits + # Filter non-zero rates for generator + nz_rates = [] + nz_generators = [] + threshold = 1e-5 + for rate, gen in zip(rates, generators): + if rate > threshold: + nz_rates.append(rate) + nz_generators.append(gen) + self._generators = nz_generators + self._rates = np.array(nz_rates, dtype=float) + + # Parameters to be initialized from generators and rates + self._generator_mats = {} + self._noise_strengths = {} + self._sampling_mats = {} + + # RNG for CTMP sampling + self._rng = None + self.seed(seed) + + def expectation_value(self, + counts: Dict, + diagonal: Optional[np.ndarray] = None, + qubits: Optional[List[int]] = None, + clbits: Optional[List[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: + counts: counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the 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. + + Returns: + (float, float): the expectation value and 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 = list(range(self._num_qubits)) + + # Convert counts to probs + probs, shots = counts_probability_vector( + counts, clbits=clbits, qubits=qubits, + num_qubits=len(qubits), return_shots=True) + + # Ensure diagonal is a numpy vector so we can use fancy indexing + if diagonal is None: + diagonal = self._z_diagonal(2**len(qubits)) + diagonal = np.asarray(diagonal, dtype=probs.dtype) + + # Get arrays from CSC sparse matrix format + # TODO: To use qubits kwarg we should compute the B-matrix and gamma + # for the specified qubit subsystem + gamma = self.noise_strength(qubits) + bmat = self._sampling_matrix(qubits) + values = bmat.data + indices = np.asarray(bmat.indices, dtype=np.int) + indptrs = np.asarray(bmat.indptr, dtype=np.int) + + # Set a minimum number of CTMP samples + shots = sum(counts.values()) + min_delta = 0.05 + shots_delta = max(4 / (min_delta**2), shots) + num_samples = int(np.ceil(shots_delta * np.exp(2 * gamma))) + + # Break total number of samples up into steps of a max number + # of samples + expval = 0 + batch_size = 50000 + samples_set = (num_samples // batch_size) * [batch_size] + [ + num_samples % batch_size] + for sample_shots in samples_set: + # Apply sampling + samples, sample_signs = self._ctmp_inverse( + sample_shots, probs, gamma, values, indices, indptrs, self._rng) + + # Compute expectation value + expval += diagonal[samples[sample_signs == 0]].sum() + expval -= diagonal[samples[sample_signs == 1]].sum() + + expval = (np.exp(2 * gamma) / num_samples) * expval + + # TODO: calculate exact standard deviation + # For now we return the upper bound: stddev <= Gamma / sqrt(shots) + # using Gamma ~ exp(2 * gamma) + stddev = np.exp(2 * gamma) / np.sqrt(shots) + return expval, stddev + + def generator_matrix(self, qubits: List[int] = None) -> sps.coo_matrix: + r"""Return the generator matrix on the specified qubits. + + The generator matrix :math:`G` is given by :math:`\sum_i r_i G_i` + where the sum is taken over all :math:`G_i` acting on the specified + qubits subset. + + Args: + qubits: Optional, qubit subset for the generators. + + Returns: + sps.coo_matrix: the generator matrix :math:`G`. + """ + if qubits is None: + qubits = tuple(range(self._num_qubits)) + else: + qubits = tuple(sorted(qubits)) + + if qubits not in self._generator_mats: + # Construct G from subset generators + qubits_set = set(qubits) + g_mat = sps.coo_matrix(2 * (2**len(qubits),), dtype=np.float) + for gen, rate in zip(self._generators, self._rates): + if qubits_set.issuperset(gen[2]): + # Keep generator + g_mat += rate * self._generator_to_coo_matrix(gen, qubits) + + # Add generator matrix to cache + self._generator_mats[qubits] = sps.coo_matrix(g_mat) + return self._generator_mats[qubits] + + def mitigation_matrix(self, qubits: List[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. + + Returns: + np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`. + """ + # NOTE: the matrix definition of G is somehow flipped in both row and + # columns compared to the canonical ordering for the A-matrix used + # in the Complete and Tensored methods + gmat = self.generator_matrix(qubits) + gmat = np.flip(gmat.todense()) + return la.expm(-gmat) + + 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. + """ + # NOTE: the matrix definition of G is somehow flipped in both row and + # columns compared to the canonical ordering for the A-matrix used + # in the Complete and Tensored methods + gmat = self.generator_matrix(qubits) + gmat = np.flip(gmat.todense()) + return la.expm(gmat) + + def noise_strength(self, qubits: Optional[int] = None) -> float: + """Return the noise strength :math:`gamma` on the specified qubits""" + if qubits is None: + qubits = tuple(range(self._num_qubits)) + else: + qubits = tuple(sorted(qubits)) + + if qubits not in self._noise_strengths: + # Compute gamma and cache + g_mat = self.generator_matrix(qubits) + # Check ideal case + if g_mat.row.size == 0: + gamma = 0 + else: + gamma = np.max(-g_mat.data[g_mat.row == g_mat.col]) + if gamma < 0: + raise QiskitError( + 'gamma should be non-negative, found gamma={}'.format(gamma)) + self._noise_strengths[qubits] = gamma + return self._noise_strengths[qubits] + + def seed(self, value=None): + """Set the seed for the quantum state RNG.""" + if isinstance(value, np.random.Generator): + self._rng = value + else: + self._rng = np.random.default_rng(value) + + def _compute_gamma(self, qubits=None): + """Compute gamma for N-qubit mitigation""" + if qubits is not None: + raise NotImplementedError( + "qubits kwarg is not yet implemented for CTMP method.") + gamma = self.noise_strength(qubits) + return np.exp(2 * gamma) + + @staticmethod + def _ctmp_inverse( + n_samples: int, + probs: np.ndarray, + gamma: float, + csc_data: np.ndarray, + csc_indices: np.ndarray, + csc_indptrs: np.ndarray, + rng: np.random.Generator) -> Tuple[Tuple[int], Tuple[int]]: + """Apply CTMP algorithm to input counts dictionary, return + sampled counts and associated shots. Equivalent to Algorithm 1 in + Bravyi et al. + + Args: + n_samples: Number of times to sample in CTMP algorithm. + probs: probability vector constructed from counts. + gamma: noise strength parameter + csc_data: Sparse CSC matrix data array (`csc_matrix.data`). + csc_indices: Sparse CSC matrix indices array (`csc_matrix.indices`). + csc_indptrs: Sparse CSC matrix indices array (`csc_matrix.indptrs`). + rng: RNG generator object. + + Returns: + Tuple[Tuple[int], Tuple[int]]: Resulting list of shots and associated + signs from the CTMP algorithm. + """ + alphas = rng.poisson(lam=gamma, size=n_samples) + signs = np.mod(alphas, 2) + x_vals = rng.choice(len(probs), size=n_samples, p=probs) + + # Apply CTMP sampling + r_vals = rng.random(size=alphas.sum()) + y_vals = np.zeros(x_vals.size, dtype=np.int) + _markov_chain_compiled(y_vals, x_vals, r_vals, alphas, csc_data, csc_indices, + csc_indptrs) + + return y_vals, signs + + def _sampling_matrix(self, qubits: Optional[int] = None) -> sps.csc_matrix: + """Compute the B matrix for CTMP sampling""" + if qubits is None: + qubits = tuple(range(self._num_qubits)) + else: + qubits = tuple(sorted(qubits)) + + if qubits not in self._sampling_mats: + # Compute matrix and cache + gmat = self.generator_matrix(qubits) + gamma = self.noise_strength(qubits) + bmat = sps.eye(2**len(qubits)) + if gamma != 0: + bmat = bmat + gmat / gamma + self._sampling_mats[qubits] = bmat.tocsc() + + return self._sampling_mats[qubits] + + @staticmethod + def _tensor_list(parts: List[np.ndarray]) -> np.ndarray: + """Compute sparse tensor product of all matrices in list""" + res = parts[0] + for mat in parts[1:]: + res = sps.kron(res, mat) + return res + + def _generator_to_coo_matrix(self, gen: Generator, qubits: Tuple[int]) -> sps.coo_matrix: + """Compute the matrix form of a generator. + + Generators are uniquely determined by two bitstrings, + and a list of qubits on which the bitstrings act. For instance, + the generator `|b^i> int: + """Choise a random array element from specified distribution. + + Given a list and associated probabilities for each element of the list, + return a random choice. This function is required since Numpy + random.choice cannot be compiled with Numba. + + Args: + inds (List[int]): List of indices to choose from. + probs (List[float]): List of probabilities for indices. + r_val: float: pre-generated random number in [0, 1). + + Returns: + int: Randomly sampled index from list. + """ + probs = np.cumsum(probs) + n_probs = len(probs) + for i in range(n_probs): + if r_val < probs[i]: + return inds[i] + return inds[-1] + + +@jit_fallback +def _markov_chain_compiled(y_vals: np.ndarray, x_vals: np.ndarray, + r_vals: np.ndarray, alpha_vals: np.ndarray, + csc_vals: np.ndarray, csc_indices: np.ndarray, + csc_indptrs: np.ndarray): + """Simulate the Markov process for a CSC transition matrix. + + Args: + y_vals: array to store sampled values. + x_vals: array of initial state values. + r_vals: pre-generated array of random numbers [0, 1) to use in sampling. + alpha_vals: array of Markov step values for sampling. + csc_vals: Sparse CSC matrix data array (`csc_matrix.data`). + csc_indices: Sparse CSC matrix indices array (`csc_matrix.indices`). + csc_indptrs: Sparse CSC matrix indices array (`csc_matrix.indptrs`). + """ + num_samples = y_vals.size + r_pos = 0 + for i in range(num_samples): + y = x_vals[i] + for _ in range(alpha_vals[i]): + begin_slice = csc_indptrs[y] + end_slice = csc_indptrs[y + 1] + probs = csc_vals[begin_slice:end_slice] + inds = csc_indices[begin_slice:end_slice] + y = _choice(inds, probs, r_vals[r_pos]) + r_pos += 1 + y_vals[i] = y diff --git a/qiskit/ignis/mitigation/expval/fitter.py b/qiskit/ignis/mitigation/expval/fitter.py new file mode 100644 index 000000000..ff7824fd0 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/fitter.py @@ -0,0 +1,89 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. +""" +Expectation value measurement migitator fitter. +""" + +from typing import Optional, Dict, List, Union + +from qiskit.result import Result +from qiskit.exceptions import QiskitError + +from .utils import calibration_data, assignment_matrix +from .complete_mitigator import CompleteExpvalMeasMitigator +from .tensored_mitigator import TensoredExpvalMeasMitigator +from .ctmp_mitigator import CTMPExpvalMeasMitigator +from .ctmp_fitter import fit_ctmp_meas_mitigator +from .ctmp_generator_set import Generator + + +class ExpvalMeasMitigatorFitter: + """Expectation value measurement error mitigator calibration fitter. + + See :func:`qiskit.ignis.mitigation.expval_meas_mitigator_circuits` for + additional documentation. + """ + + def __init__(self, + result: Result, + metadata: List[Dict[str, any]]): + """Fit a measurement error mitigator object from experiment data. + + Args: + result: Qiskit result object. + metadata: mitigation generator metadata. + """ + self._num_qubits = None + self._cal_data = None + self._mitigator = None + self._cal_data, self._num_qubits, self._method = calibration_data( + result, metadata) + + @property + def mitigator(self): + """Return the fitted mitigator object""" + if self._mitigator is None: + raise QiskitError("Mitigator has not been fitted. Run `fit` first.") + return self._mitigator + + def fit(self, method: Optional[str] = None, + generators: Optional[List[Generator]] = None) -> Union[ + CompleteExpvalMeasMitigator, + TensoredExpvalMeasMitigator, + CTMPExpvalMeasMitigator]: + """Fit and return the Mitigator object from the calibration data.""" + + if method is None: + method = self._method + + if method == 'complete': + # Construct A-matrix from calibration data + amat = assignment_matrix(self._cal_data, self._num_qubits) + self._mitigator = CompleteExpvalMeasMitigator(amat) + + elif method == 'tensored': + # Construct single-qubit A-matrices from calibration data + amats = [] + for qubit in range(self._num_qubits): + amat = assignment_matrix(self._cal_data, self._num_qubits, [qubit]) + amats.append(amat) + self._mitigator = TensoredExpvalMeasMitigator(amats) + + elif method in ['CTMP', 'ctmp']: + self._mitigator = fit_ctmp_meas_mitigator( + self._cal_data, self._num_qubits, generators) + else: + raise QiskitError( + "Invalid expval measurement error mitigation method {}".format(method)) + return self._mitigator diff --git a/qiskit/ignis/mitigation/expval/tensored_mitigator.py b/qiskit/ignis/mitigation/expval/tensored_mitigator.py new file mode 100644 index 000000000..e20350063 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/tensored_mitigator.py @@ -0,0 +1,195 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +Single-qubit tensor-product measurement error mitigation generator. +""" +from typing import Optional, List, Dict, Tuple +import numpy as np + +from .utils import (counts_probability_vector, _expval_with_stddev) +from .base_meas_mitigator import BaseExpvalMeasMitigator + + +class TensoredExpvalMeasMitigator(BaseExpvalMeasMitigator): + """1-qubit tensor product measurement error mitigator. + + This class can be used with the + :func:`qiskit.ignis.mitigation.expectation_value` function to apply + measurement error mitigation of local single-qubit measurement errors. + Expectation values can also be computed directly using the + :meth:`expectation_value` method. + + For measurement mitigation to be applied the mitigator should be + calibrated using the + :func:`qiskit.ignis.mitigation.expval_meas_mitigator_circuits` function + and :class:`qiskit.ignis.mitigation.ExpvalMeasMitigatorFitter` class with + the ``'tensored'`` mitigation method. + """ + + def __init__(self, amats: List[np.ndarray]): + """Initialize a TensorMeasurementMitigator + + Args: + amats: list of single-qubit readout error assignment matrices. + """ + self._num_qubits = len(amats) + self._assignment_mats = amats + 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 + + def expectation_value(self, + counts: Dict, + diagonal: Optional[np.ndarray] = None, + qubits: Optional[List[int]] = None, + clbits: Optional[List[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: + counts: counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the 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. + + Returns: + (float, float): the expectation value and 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)``. + """ + # Get expectation value on specified qubits + probs, shots = counts_probability_vector( + counts, clbits=clbits, return_shots=True) + num_qubits = int(np.log2(probs.shape[0])) + + if qubits is not None: + ainvs = self._mitigation_mats[list(qubits)] + else: + ainvs = self._mitigation_mats + + # Get operator coeffs + if diagonal is None: + diagonal = self._z_diagonal(2 ** num_qubits) + # 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() + + return _expval_with_stddev(coeffs, probs, shots) + + def mitigation_matrix(self, qubits: List[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. + + Returns: + np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`. + """ + if qubits is None: + qubits = list(range(self._num_qubits)) + if isinstance(qubits, int): + qubits = [qubits] + mat = self._mitigation_mats[qubits[0]] + for i in qubits[1:]: + mat = np.kron(self._mitigation_mats[qubits[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 = list(range(self._num_qubits)) + if isinstance(qubits, int): + qubits = [qubits] + mat = self._assignment_mats[qubits[0]] + for i in qubits[1:]: + mat = np.kron(self._assignment_mats[qubits[i]], mat) + return mat + + def assignment_fidelity(self, qubits: Optional[List[int]] = None) -> float: + r"""Return the measurement assignment fidelity on the specified qubits. + + The assignment fidelity on N-qubits is defined as + :math:`\sum_{x\in\{0, 1\}^n} P(x|x) / 2^n`, where + :math:`P(x|x) = \rangle x|A|x\langle`, and :math:`A` is the + :meth:`assignment_matrix`. + + Args: + qubits: Optional, qubits being measured for operator expval. + + Returns: + float: the assignment fidelity. + """ + if qubits is None: + qubits = list(range(self._num_qubits)) + if isinstance(qubits, int): + qubits = [qubits] + fid = 1.0 + for i in qubits: + fid *= np.mean(self._assignment_mats[i].diagonal()) + return fid + + def _compute_gamma(self, qubits=None): + """Compute gamma for N-qubit mitigation""" + if qubits is None: + gammas = self._gammas + else: + gammas = self._gammas[list(qubits)] + return np.product(gammas) diff --git a/qiskit/ignis/mitigation/expval/utils.py b/qiskit/ignis/mitigation/expval/utils.py new file mode 100644 index 000000000..bb9bb42a0 --- /dev/null +++ b/qiskit/ignis/mitigation/expval/utils.py @@ -0,0 +1,292 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. +""" +Measurement error mitigation utility functions. +""" + +import logging +from functools import partial +from typing import Optional, List, Dict, Tuple +import numpy as np + +from qiskit.exceptions import QiskitError +from qiskit.result import Counts, Result +from qiskit.ignis.verification.tomography import marginal_counts, combine_counts +from qiskit.ignis.numba import jit_fallback + +logger = logging.getLogger(__name__) + + +def expectation_value(counts: Counts, + diagonal: Optional[np.ndarray] = None, + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, + meas_mitigator: Optional = None, + ) -> Tuple[float, float]: + r"""Compute the expectation value of a diagonal operator from counts. + + This computes the estimator of + :math:`\langle O \rangle = \mbox{Tr}[\rho. O]`, optionally with measurement + error mitigation, of a diagonal observable + :math:`O = \sum_{x\in\{0, 1\}^n} O(x)|x\rangle\!\langle x|`. + + Args: + counts: counts object + diagonal: Optional, the vector of diagonal values for summing the + expectation value. If ``None`` the 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. + meas_mitigator: Optional, a measurement mitigator to apply mitigation. + + Returns: + (float, float): the expectation value and 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)``. + + For calibrating a expval measurement error mitigator for the + ``meas_mitigator`` kwarg see + :func:`qiskit.ignis.mitigation.expval_meas_mitigator_circuits` and + :class:`qiskit.ignis.mitigation.ExpvalMeasMitigatorFitter`. + """ + if meas_mitigator is not None: + # Use mitigator expectation value method + return meas_mitigator.expectation_value( + counts, diagonal=diagonal, clbits=clbits, + qubits=qubits) + + # Marginalize counts + if clbits is not None: + counts = marginal_counts(counts, meas_qubits=clbits) + + # Get counts shots and probabilities + probs = np.array(list(counts.values())) + shots = probs.sum() + probs = probs / shots + + # Get diagonal operator coefficients + if diagonal is None: + coeffs = np.array([(-1) ** (key.count('1') % 2) + for key in counts.keys()], dtype=probs.dtype) + else: + diagonal = np.asarray(diagonal) + keys = [int(key, 2) for key in counts.keys()] + coeffs = np.asarray(diagonal[keys], dtype=probs.dtype) + + return _expval_with_stddev(coeffs, probs, shots) + + +def counts_probability_vector( + counts: Counts, + qubits: Optional[List[int]] = None, + clbits: Optional[List[int]] = None, + num_qubits: Optional[int] = None, + return_shots: Optional[bool] = False) -> np.ndarray: + """Compute mitigated expectation value. + + Args: + counts: counts object + qubits: qubits the count bitstrings correspond to. + clbits: Optional, marginalize counts to just these bits. + num_qubits: the total number of qubits. + return_shots: return the number of shots. + + Raises: + QiskitError: if qubit and clbit kwargs are not valid. + + Returns: + np.ndarray: a probability vector for all count outcomes. + """ + # Marginalize counts + if clbits is not None: + counts = marginal_counts(counts, meas_qubits=clbits) + + # Get total number of qubits + if num_qubits is None: + num_qubits = len(next(iter(counts))) + + # Get 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 + + # Remap 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) + if return_shots: + return vec, shots + return vec + + +def calibration_data(result: Result, + metadata: List[Dict[str, any]]) -> Tuple[ + Dict[int, Dict[int, int]], int]: + """Return FullMeasureErrorMitigator from result data. + + Args: + result: Qiskit result object. + metadata: mitigation generator metadata. + + Returns: + Calibration data dictionary {label: Counts} and number of qubits. + """ + # Filter mitigation calibration data + cal_data = {} + num_qubits = None + method = None + for i, meta in enumerate(metadata): + if meta.get('experiment') == 'meas_mit': + if num_qubits is None: + num_qubits = len(meta['cal']) + if method is None: + method = meta.get('method', None) + key = int(meta['cal'], 2) + counts = result.get_counts(i).int_outcomes() + if key not in cal_data: + cal_data[key] = counts + else: + cal_data[key] = combine_counts(cal_data[key], counts) + return cal_data, num_qubits, method + + +def assignment_matrix(cal_data: Dict[int, Dict[int, int]], + num_qubits: int, + qubits: Optional[List[int]] = None) -> np.array: + """Computes the assignment matrix for specified qubits. + + Args: + cal_data: calibration dataset. + num_qubits: the number of qubits for the calibation dataset. + qubits: Optional, the qubit subset to construct the matrix on. + + Returns: + np.ndarray: the constructed A-matrix. + + Raises: + QiskitError: if the calibration data is not sufficient for + reconstruction on the specified qubits. + """ + # If qubits is None construct full A-matrix from calibration data + # Otherwise we compute the local A-matrix on specified + # qubits subset. This involves filtering the cal data + # on no-errors occuring on qubits outside the subset + + if qubits is not None: + qubits = np.asarray(qubits) + dim = 1 << qubits.size + mask = _amat_mask(qubits, num_qubits) + accum_func = partial(_amat_accum_local, qubits, mask) + + else: + qubits = np.array(range(num_qubits)) + dim = 1 << num_qubits + accum_func = _amat_accum_full + + amat = np.zeros([dim, dim], dtype=float) + for cal, counts in cal_data.items(): + counts_keys = np.array(list(counts.keys())) + counts_values = np.array(list(counts.values())) + accum_func(amat, cal, counts_keys, counts_values) + + renorm = amat.sum(axis=0, keepdims=True) + if np.any(renorm == 0): + raise QiskitError( + 'Insufficient calibration data to fit assignment matrix' + ' on qubits {}'.format(qubits.tolist())) + return amat / renorm + + +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) + stddev = np.sqrt(variance) if variance > 0 else 0.0 + return expval, stddev + + +@jit_fallback +def _amat_mask(qubits, num_qubits): + """Compute bit-mask for other other non-specified qubits.""" + mask = 0 + for i in range(num_qubits): + if not np.any(qubits == i): + mask += 1 << i + return mask + + +@jit_fallback +def _amat_index(i, qubits): + """Compute local index for specified qubits and full index value.""" + masks = 1 << qubits + shifts = np.arange(qubits.size - 1, -1, -1) + val = ((i & masks) >> qubits) << shifts + return np.sum(val) + + +@jit_fallback +def _amat_accum_local(qubits, mask, mat, cal, counts_keys, counts_values): + """Accumulate calibration data on the specified matrix""" + x_out = cal & mask + x_ind = _amat_index(cal, qubits) + for i in range(counts_keys.size): + y_out = counts_keys[i] & mask + if x_out == y_out: + y_ind = _amat_index(counts_keys[i], qubits) + mat[y_ind, x_ind] += counts_values[i] + + +@jit_fallback +def _amat_accum_full(mat, cal, counts_keys, counts_values): + """Accumulate calibration data on the specified matrix""" + for i in range(counts_keys.size): + mat[counts_keys[i], cal] += counts_values[i] diff --git a/qiskit/ignis/mitigation/measurement/circuits.py b/qiskit/ignis/mitigation/measurement/circuits.py index dabe8911b..d05a9ee74 100644 --- a/qiskit/ignis/mitigation/measurement/circuits.py +++ b/qiskit/ignis/mitigation/measurement/circuits.py @@ -19,7 +19,7 @@ from typing import List, Tuple, Union from qiskit import QuantumRegister, ClassicalRegister, \ QuantumCircuit, QiskitError -from ...verification.tomography import count_keys +from qiskit.ignis.verification.tomography import count_keys def complete_meas_cal(qubit_list: List[int] = None, diff --git a/qiskit/ignis/mitigation/measurement/filters.py b/qiskit/ignis/mitigation/measurement/filters.py index 7f5abf8e1..4366bdd75 100644 --- a/qiskit/ignis/mitigation/measurement/filters.py +++ b/qiskit/ignis/mitigation/measurement/filters.py @@ -26,7 +26,7 @@ import qiskit from qiskit import QiskitError from qiskit.tools import parallel_map -from ...verification.tomography import count_keys +from qiskit.ignis.verification.tomography import count_keys class MeasurementFilter(): diff --git a/qiskit/ignis/mitigation/measurement/fitters.py b/qiskit/ignis/mitigation/measurement/fitters.py index 5371a19ba..cef15a4da 100644 --- a/qiskit/ignis/mitigation/measurement/fitters.py +++ b/qiskit/ignis/mitigation/measurement/fitters.py @@ -24,8 +24,9 @@ import numpy as np from qiskit import QiskitError from qiskit.result import Result +from qiskit.ignis.verification.tomography import count_keys from .filters import MeasurementFilter, TensoredFilter -from ...verification.tomography import count_keys + try: from matplotlib import pyplot as plt diff --git a/qiskit/ignis/numba.py b/qiskit/ignis/numba.py new file mode 100644 index 000000000..88dfa426a --- /dev/null +++ b/qiskit/ignis/numba.py @@ -0,0 +1,35 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2020. +# +# 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. + +"""Optional support for Numba just-in-time compilation.""" + +import logging +logger = logging.getLogger(__name__) + +try: + import numba + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + logger.info('Numba not installed, certain functions have improved ' + 'performance if Number is installed: ' + 'https://pypi.org/project/numba/') + + +def jit_fallback(func): + """Decorator to try to apply numba JIT compilation.""" + if _HAS_NUMBA: + return numba.jit(nopython=True)(func) + else: + return func diff --git a/qiskit/ignis/verification/tomography/data.py b/qiskit/ignis/verification/tomography/data.py index a8795f2f9..93e5ccb34 100644 --- a/qiskit/ignis/verification/tomography/data.py +++ b/qiskit/ignis/verification/tomography/data.py @@ -73,7 +73,7 @@ def marginal_counts(counts: Dict[str, int], # Generate bitstring keys for measured qubits meas_keys = count_keys(len(qubits)) - # Get regex match strings for suming outcomes of other qubits + # Get regex match strings for summing outcomes of other qubits rgx = [] for key in meas_keys: def helper(x, y): diff --git a/releasenotes/notes/expval-meas-mit-0f3d97311ef80aa9.yaml b/releasenotes/notes/expval-meas-mit-0f3d97311ef80aa9.yaml new file mode 100644 index 000000000..a8b02db0a --- /dev/null +++ b/releasenotes/notes/expval-meas-mit-0f3d97311ef80aa9.yaml @@ -0,0 +1,79 @@ +--- +features: + - | + Adds expectation value measurement error mitigation to the mitigation module. + This supports using *complete* N-qubit assignment matrix, single-qubit + *tensored* assignment matrix, or *continuous time Markov process (CTMP)* [1] + measurement error mitigation when computing expectation values of diagonal + operators from counts dictionaries. Expectation values are computed using + the using the :func:`qiskit.ignis.mitigation.expectation_value` function. + + Calibration circuits for calibrating a measurement error mitigator are + generated using the :func:`qiskit.ignis.mitigation.expval_meas_mitigator_circuits` + function, and the result fitted using the + :class:`qiskit.ignis.mitigation.ExpvalMeasMitigatorFitter` class. The + fitter returns a mitigator object can the be supplied as an argument to the + :func:`~qiskit.ignis.mitigation.expectation_value` function to apply mitigation. + + [1] S Bravyi, S Sheldon, A Kandala, DC Mckay, JM Gambetta, + *Mitigating measurement errors in multi-qubit experiments*, + arXiv:2006.14044 [quant-ph]. + + Example: + + The following example shows calibrating a 5-qubit expectation value + measurement error mitigator using the ``'tensored'`` method. + + .. jupyter-execute:: + + from qiskit import execute + from qiskit.test.mock import FakeVigo + import qiskit.ignis.mitigation as mit + + backend = FakeVigo() + num_qubits = backend.configuration().num_qubits + + # Generate calibration circuits + circuits, metadata = mit.expval_meas_mitigator_circuits( + num_qubits, method='tensored') + result = execute(circuits, backend, shots=8192).result() + + # Fit mitigator + mitigator = mit.ExpvalMeasMitigatorFitter(result, metadata).fit() + + # Plot fitted N-qubit assignment matrix + mitigator.plot_assignment_matrix() + + The following shows how to use the above mitigator to apply measurement + error mitigation to expectation value computations + + .. jupyter-execute:: + + from qiskit import QuantumCircuit + + # Test Circuit with expectation value -1. + qc = QuantumCircuit(num_qubits) + qc.x(range(num_qubits)) + qc.measure_all() + + # Execute + shots = 8192 + seed_simulator = 1999 + result = execute(qc, backend, shots=8192, seed_simulator=1999).result() + counts = result.get_counts(0) + + # Expectation value of Z^N without mitigation + expval_nomit, error_nomit = mit.expectation_value(counts) + print('Expval (no mitigation): {:.2f} \u00B1 {:.2f}'.format( + expval_nomit, error_nomit)) + + # Expectation value of Z^N with mitigation + expval_mit, error_mit = mit.expectation_value(counts, + meas_mitigator=mitigator) + print('Expval (with mitigation): {:.2f} \u00B1 {:.2f}'.format( + expval_mit, error_mit)) + - | + Adds Numba as an optional dependency. Numba is used to significantly increase + the performance of the :class:`qiskit.ignis.mitigation.CTMPExpvalMeasMitigator` + class used for expectation value measurement error mitigation with the CTMP + method. diff --git a/setup.py b/setup.py index 6fa4209d9..a33a92c77 100644 --- a/setup.py +++ b/setup.py @@ -79,6 +79,7 @@ extras_require={ 'visualization': ['matplotlib>=2.1'], 'cvx': ['cvxpy>=1.0.15'], + 'jit': ['numba'], }, install_requires=requirements, include_package_data=True, diff --git a/test/measurement/discriminator/test_iq_discriminator.py b/test/measurement/discriminator/test_iq_discriminator.py index a107a031b..ef4ab87e5 100644 --- a/test/measurement/discriminator/test_iq_discriminator.py +++ b/test/measurement/discriminator/test_iq_discriminator.py @@ -30,7 +30,7 @@ from qiskit.ignis.measurement.discriminator.filters import DiscriminationFilter from qiskit.ignis.measurement.discriminator.iq_discriminators import \ LinearIQDiscriminator, SklearnIQDiscriminator -from qiskit.ignis.mitigation.measurement import circuits +from qiskit.ignis.mitigation.measurement import tensored_meas_cal from qiskit.result.models import ExperimentResultData @@ -48,7 +48,7 @@ def setUp(self): self.shots = 52 self.qubits = [0, 1] - meas_cal, _ = circuits.tensored_meas_cal([[0], [1]]) + meas_cal, _ = tensored_meas_cal([[0], [1]]) backend = Aer.get_backend('qasm_simulator') job = qiskit.execute(meas_cal, backend=backend, shots=self.shots, diff --git a/test/mitigation/__init__.py b/test/mitigation/__init__.py new file mode 100644 index 000000000..d2a3dbeec --- /dev/null +++ b/test/mitigation/__init__.py @@ -0,0 +1,13 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. diff --git a/test/mitigation/test_assignment_matrices.py b/test/mitigation/test_assignment_matrices.py new file mode 100644 index 000000000..9066af703 --- /dev/null +++ b/test/mitigation/test_assignment_matrices.py @@ -0,0 +1,127 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +Test assignment matrices +""" + +import unittest +from typing import List +from functools import lru_cache +from itertools import product + +from ddt import ddt, unpack, data +import numpy as np + +from qiskit import QuantumCircuit +from qiskit.ignis.mitigation import ( + expval_meas_mitigator_circuits, + ExpvalMeasMitigatorFitter +) + +from .test_mitigators import NoisySimulationTest + + +_QUBIT_SUBSETS = [(1, 0, 2), (3, 0), (2,), (0, 3, 1, 2)] +_PAIRINGS = [ + *product(['complete', 'tensored'], _QUBIT_SUBSETS, [True, False]), + ('ctmp', None, True), + ('ctmp', None, False), +] + + +@ddt +class TestMatrices(NoisySimulationTest): + """Test that assignment matrices are computed correctly (close to identity) + and produce correct expectation values. + """ + + @lru_cache(6) + def get_mitigator(self, method: str, noise_model: bool): + """Return the mitigator for the given parameters. + """ + circs, meta = expval_meas_mitigator_circuits(self.num_qubits, method=method) + cal_res = self.execute_circs(circs, noise_model=self.noise_model if noise_model else None) + mitigator = ExpvalMeasMitigatorFitter(cal_res, meta).fit(method=method) + return mitigator + + @data(*_PAIRINGS) + @unpack + def test_matrices(self, method: str, qubits: List[int], noise: bool): + """Compute and compare matrices with given options. + """ + if qubits is not None: + num_qubits = len(qubits) + else: + num_qubits = self.num_qubits + mitigator = self.get_mitigator(method, noise_model=noise) + assignment_matrix = mitigator.assignment_matrix(qubits) + mitigation_matrix = mitigator.mitigation_matrix(qubits) + + np.testing.assert_array_almost_equal( + assignment_matrix, + np.eye(2**num_qubits, 2**num_qubits), + decimal=1 if noise else 6 + ) + + np.testing.assert_array_almost_equal( + mitigation_matrix, + np.eye(2**num_qubits, 2**num_qubits), + decimal=1 if noise else 6 + ) + + @data( + *product( + ['complete', 'tensored'], + ['1011', '1010', '1111'], + [[0, 1, 2], [1], [1, 0, 3, 2], None] + ), + *product( + ['ctmp'], + ['1011', '1010', '1111'], + [None] + ) + ) + @unpack + def test_parity_exp_vals_partial(self, method: str, bitstr: str, qubits: List[int]): + """Compute expectation value of parity operators and + compare with exact result. + """ + mitigator = self.get_mitigator(method, True) + + if qubits is None: + exp_val_exact = (-1)**(bitstr.count('1')) + else: + exp_val_exact = 1 + for q in qubits: + exp_val_exact *= (-1)**int(bitstr[::-1][q]) + + qc = QuantumCircuit(self.num_qubits, self.num_qubits) + for i, b in enumerate(bitstr[::-1]): + if b == '1': + qc.x(i) + qc.measure(range(4), range(4)) + + counts = self.execute_circs(qc, noise_model=self.noise_model).get_counts(0) + + exp_val_mit, _ = mitigator.expectation_value( + counts=counts, + clbits=qubits, + qubits=qubits + ) + + self.assertLessEqual(abs(exp_val_exact - exp_val_mit), 0.1) + + +if __name__ == '__main__': + unittest.main() diff --git a/test/mitigation/test_markov.py b/test/mitigation/test_markov.py new file mode 100644 index 000000000..fc6a4324b --- /dev/null +++ b/test/mitigation/test_markov.py @@ -0,0 +1,115 @@ + +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +CTMP markov sampling tests +""" +import unittest + +import numpy as np +import scipy as sp +from qiskit.ignis.mitigation import CTMPExpvalMeasMitigator +from qiskit.ignis.mitigation.expval.ctmp_mitigator import _markov_chain_compiled + + +def statistical_test(num_tests: int, fraction_passes: float): + """Wrapper for statistical test that should be run multiple times and pass + at least a certain fraction of times. + """ + # pylint: disable=broad-except + def stat_test(func): + def wrapper_func(*args, **kwargs): + + num_failures = 0 + num_passes = 0 + for _ in range(num_tests): + try: + func(*args, **kwargs) + num_passes += 1 + except Exception: + num_failures += 1 + if num_passes / num_tests < fraction_passes: + raise ValueError('Passed {} out of {} trials, needed {}%'.format( + num_passes, + num_tests, + 100 * fraction_passes + )) + + return wrapper_func + + return stat_test + + +class TestMarkov(unittest.TestCase): + """Test the (un)compiled markov process simulator + """ + + def setUp(self): + self.r_dict = { + # (final, start, qubits) + ('0', '1', (0,)): 1e-3, + ('1', '0', (0,)): 1e-1, + ('0', '1', (1,)): 1e-3, + ('1', '0', (1,)): 1e-1, + + ('00', '11', (0, 1)): 1e-3, + ('11', '00', (0, 1)): 1e-1, + ('01', '10', (0, 1)): 1e-3, + ('10', '01', (0, 1)): 1e-3 + } + + mitigator = CTMPExpvalMeasMitigator( + generators=self.r_dict.keys(), + rates=self.r_dict.values(), + num_qubits=2 + ) + self.trans_mat = sp.linalg.expm(mitigator.generator_matrix()).tocsc() + self.num_steps = 100 + + def markov_chain_int(self, x): + """Convert input for markov_chain_compiled function""" + y = np.array([0]) + x = np.array([x]) + alpha = np.array([self.num_steps]) + rng = np.random.default_rng(seed=x) + r_vals = rng.random(size=alpha.sum()) + values = self.trans_mat.data + indices = np.asarray(self.trans_mat.indices, dtype=np.int) + indptrs = np.asarray(self.trans_mat.indptr, dtype=np.int) + _markov_chain_compiled(y, x, r_vals, alpha, values, indices, indptrs) + return y[0] + + @statistical_test(50, 0.7) + def test_markov_chain_int_0(self): + """Test markov process starting at specific state""" + self.assertEqual(self.markov_chain_int(0), 3) + + @statistical_test(50, 0.7) + def test_markov_chain_int_1(self): + """Test markov process starting at specific state""" + self.assertEqual(self.markov_chain_int(1), 3) + + @statistical_test(50, 0.7) + def test_markov_chain_int_2(self): + """Test markov process starting at specific state""" + self.assertEqual(self.markov_chain_int(2), 3) + + @statistical_test(50, 0.7) + def test_markov_chain_int_3(self): + """Test markov process starting at specific state""" + self.assertEqual(self.markov_chain_int(3), 3) + + +if __name__ == '__main__': + unittest.main() diff --git a/test/mitigation/test_mitigators.py b/test/mitigation/test_mitigators.py new file mode 100644 index 000000000..b8abd5b27 --- /dev/null +++ b/test/mitigation/test_mitigators.py @@ -0,0 +1,164 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2019, 2020. +# +# 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. +""" +Test mitigators +""" + +from itertools import combinations, product, chain, permutations +from typing import List, Tuple +import unittest + +from ddt import ddt, unpack, data + +from qiskit import execute, QuantumCircuit +from qiskit.result import Result +from qiskit.providers.aer import QasmSimulator, noise +from qiskit.ignis.mitigation import ( + expval_meas_mitigator_circuits, + ExpvalMeasMitigatorFitter, + expectation_value +) + + +class NoisySimulationTest(unittest.TestCase): + """Base class that contains methods and attributes + for doing tests of readout error noise with flexible + readout errors. + """ + + sim = QasmSimulator() + + # Example max qubit number + num_qubits = 4 + + # Create readout errors + readout_errors = [] + for i in range(num_qubits): + p_error1 = (i + 1) * 0.002 + p_error0 = 2 * p_error1 + ro_error = noise.ReadoutError([[1 - p_error0, p_error0], [p_error1, 1 - p_error1]]) + readout_errors.append(ro_error) + # TODO: Needs 2q errors? + + # Readout Error only + noise_model = noise.NoiseModel() + for i in range(num_qubits): + noise_model.add_readout_error(readout_errors[i], [i]) + seed_simulator = 100 + + shots = 10000 + + tolerance = 0.05 + + def execute_circs(self, qc_list: List[QuantumCircuit], + noise_model=None) -> Result: + """Run circuits with the readout noise defined in this class + """ + return execute( + qc_list, + backend=self.sim, + seed_simulator=self.seed_simulator, + shots=self.shots, + noise_model=None if noise_model is None else self.noise_model, + backend_options={'method': 'density_matrix'} + ).result() + + +@ddt +class TestExpVals(NoisySimulationTest): + """Test the expectation values of all the ExpvalMeasMitigator* classes + and compare against the exact results. + """ + + def setUp(self): + qc = QuantumCircuit(self.num_qubits) + qc.h(0) + for i in range(1, self.num_qubits): + qc.cx(i - 1, i) + qc.measure_all() + result_targ = self.execute_circs(qc) + result_noise = self.execute_circs(qc, noise_model=self.noise_model) + self.counts_ideal = result_targ.get_counts(0) + self.counts_noise = result_noise.get_counts(0) + + @data(*product(['complete', 'tensored', 'CTMP'], + [None, [(-1) ** i for i in range(16)], + [0.1 * i for i in range(16)]])) + @unpack + def test_expval_mitigator(self, method, diagonal): + """Test ExpvalMeasMitigator""" + circs, meta = expval_meas_mitigator_circuits(self.num_qubits, method=method) + result_cal = self.execute_circs(circs, noise_model=self.noise_model) + mitigator = ExpvalMeasMitigatorFitter(result_cal, meta).fit() + target, _ = expectation_value(self.counts_ideal, diagonal=diagonal) + expval, _ = expectation_value(self.counts_noise, + diagonal=diagonal, + meas_mitigator=mitigator) + self.assertLess(abs(target - expval), self.tolerance) + + +# https://docs.python.org/3/library/itertools.html#recipes +def powerset(iterable): + """powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)""" + output_set = list(iterable) + return chain.from_iterable(combinations(output_set, r) for r in range(len(output_set)+1)) + + +@ddt +class TestPartialExpVals(NoisySimulationTest): + """Test taking expectation values on only subsets of qubits and permutations + of the original set of qubits. + """ + + @data(*product( + [ + item for sublist in + [*map(lambda x: list(permutations(x)), powerset(range(3)))] for item in sublist], + ['complete', 'tensored'] + )) + @unpack + def test_partial_expval(self, qubits: Tuple[int], method: str): + """Test taking the expectation value only on the set of `qubits` + with `method` for mitigation. + """ + + if len(qubits) == 0: + return None + + qc = QuantumCircuit(self.num_qubits, len(qubits)) + qc.h(qubits[0]) + for i, j in combinations(qubits, r=2): + qc.cx(i, j) + qc.measure(qubits, range(len(qubits))) + + result_targ = self.execute_circs(qc) + result_nois = self.execute_circs(qc, noise_model=self.noise_model) + + counts_targ = result_targ.get_counts(0) + counts_nois = result_nois.get_counts(0) + + exp_targ, _ = expectation_value(counts_targ) + + circs, meta = expval_meas_mitigator_circuits(self.num_qubits, method=method) + result_cal = self.execute_circs(circs, noise_model=self.noise_model) + mitigator = ExpvalMeasMitigatorFitter(result_cal, meta).fit(method=method) + + expval, _ = mitigator.expectation_value(counts_nois, qubits=qubits) + + self.assertLess(abs(exp_targ - expval), self.tolerance) + return None + + +if __name__ == '__main__': + unittest.main()