# 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.
"""
Quantum process tomography analysis
"""
from typing import List, Union, Callable
from collections import defaultdict
import numpy as np
import scipy.linalg as la
from uncertainties import ufloat
from qiskit.quantum_info import DensityMatrix, Choi, Operator
from qiskit.quantum_info.operators.base_operator import BaseOperator
from qiskit.quantum_info.operators.channel.quantum_channel import QuantumChannel
from qiskit_experiments.exceptions import AnalysisError
from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options, numpy_version
from .fitters import (
tomography_fitter_data,
postprocess_fitter,
linear_inversion,
cvxpy_linear_lstsq,
cvxpy_gaussian_lstsq,
)
[docs]
class TomographyAnalysis(BaseAnalysis):
"""Base analysis for state and process tomography experiments."""
_builtin_fitters = {
"linear_inversion": linear_inversion,
"cvxpy_linear_lstsq": cvxpy_linear_lstsq,
"cvxpy_gaussian_lstsq": cvxpy_gaussian_lstsq,
}
_cvxpy_fitters = (
cvxpy_linear_lstsq,
cvxpy_gaussian_lstsq,
)
@classmethod
def _default_options(cls) -> Options:
"""Default analysis options
Analysis Options:
measurement_basis
(:class:`~qiskit_experiments.library.tomography.basis.MeasurementBasis`):
The measurement
:class:`~qiskit_experiments.library.tomography.basis.MeasurementBasis`
to use for tomographic reconstruction when running a
:class:`~qiskit_experiments.library.tomography.StateTomography` or
:class:`~qiskit_experiments.library.tomography.ProcessTomography`.
preparation_basis
(:class:`~qiskit_experiments.library.tomography.basis.PreparationBasis`):
The preparation
:class:`~qiskit_experiments.library.tomography.basis.PreparationBasis`
to use for tomographic reconstruction for
:class:`~qiskit_experiments.library.tomography.ProcessTomography`.
fitter (str or Callable): The fitter function to use for reconstruction.
This can be a string to select one of the built-in fitters, or a callable to
supply a custom fitter function. See the `Fitter Functions` section for
additional information.
fitter_options (dict): Any addition kwarg options to be supplied to the fitter
function. For documentation of available kwargs refer to the fitter function
documentation.
rescale_positive (bool): If True rescale the state returned by the fitter
to be positive-semidefinite. See the `PSD Rescaling` section for
additional information (Default: True).
rescale_trace (bool): If True rescale the state returned by the fitter
have either trace 1 for :class:`~qiskit.quantum_info.DensityMatrix`,
or trace dim for :class:`~qiskit.quantum_info.Choi` matrices (Default: True).
measurement_qubits (Sequence[int]): Optional, the physical qubits with tomographic
measurements. If not specified will be set to ``[0, ..., N-1]`` for N-qubit
tomographic measurements.
preparation_qubits (Sequence[int]): Optional, the physical qubits with tomographic
preparations. If not specified will be set to ``[0, ..., N-1]`` for N-qubit
tomographic preparations.
target (Any): Optional, target object for fidelity comparison of the fit
(Default: None).
target_bootstrap_samples (int): Optional, number of outcome re-samples to draw
from measurement data for each basis for computing a bootstrapped standard
error of fidelity with the target state. If 0 no bootstrapping will be
performed and the target fidelity will not include a standard error
(Default: 0).
target_bootstrap_seed (int | None | Generator): Optional, RNG seed or
Generator to use for bootstrapping data for boostrapped fidelity
standard error calculation (Default: None).
conditional_circuit_clbits (list[int]): Optional, the clbit indices in the
source circuit to be conditioned on when reconstructing the state.
Enabling this will return a list of reconstructed state components
conditional on the values of these clbit values. The integer value of the
conditioning clbits is stored in state analysis result extra field
`"conditional_circuit_outcome"`.
conditional_measurement_indices (list[int]): Optional, indices of tomography
measurement qubits to used for conditional state reconstruction. Enabling
this will return a list of reconstructed state components conditioned on
the remaining tomographic bases conditional on the basis index, and outcome
value for these measurements. The conditional measurement basis index and
integer value of the measurement outcome is stored in state analysis result
extra fields `"conditional_measurement_index"` and
`"conditional_measurement_outcome"` respectively.
conditional_preparation_indices (list[int]): Optional, indices of tomography
preparation qubits to used for conditional state reconstruction. Enabling
this will return a list of reconstructed channel components conditioned on
the remaining tomographic bases conditional on the basis index. The
conditional preparation basis index is stored in state analysis result
extra fields `"conditional_preparation_index"`.
extra (Dict[str, Any]): Extra metadata dictionary attached to analysis results.
"""
options = super()._default_options()
options.measurement_basis = None
options.preparation_basis = None
options.fitter = "linear_inversion"
options.fitter_options = {}
options.rescale_positive = True
options.rescale_trace = True
options.measurement_qubits = None
options.preparation_qubits = None
options.target = None
options.target_bootstrap_samples = 0
options.target_bootstrap_seed = None
options.conditional_circuit_clbits = None
options.conditional_measurement_indices = None
options.conditional_preparation_indices = None
options.extra = {}
return options
@classmethod
def _get_fitter(cls, fitter: Union[str, Callable]) -> Callable:
"""Return fitter function for named builtin fitters"""
if fitter is None:
raise AnalysisError("No tomography fitter given")
if not isinstance(fitter, str):
return fitter
if fitter in cls._builtin_fitters:
return cls._builtin_fitters[fitter]
raise AnalysisError(f"Unrecognized tomography fitter {fitter}")
def _run_analysis(self, experiment_data):
# Get option values
meas_basis = self.options.measurement_basis
meas_qubits = self.options.measurement_qubits
if meas_basis and meas_qubits is None:
meas_qubits = experiment_data.metadata.get("m_qubits")
prep_basis = self.options.preparation_basis
prep_qubits = self.options.preparation_qubits
if prep_basis and prep_qubits is None:
prep_qubits = experiment_data.metadata.get("p_qubits")
cond_meas_indices = self.options.conditional_measurement_indices
if cond_meas_indices is True:
cond_meas_indices = list(range(len(meas_qubits)))
cond_prep_indices = self.options.conditional_preparation_indices
if cond_prep_indices is True:
cond_prep_indices = list(range(len(prep_qubits)))
# Generate tomography fitter data
outcome_shape = None
if meas_basis and meas_qubits:
outcome_shape = meas_basis.outcome_shape(meas_qubits)
outcome_data, shot_data, meas_data, prep_data = tomography_fitter_data(
experiment_data.data(),
outcome_shape=outcome_shape,
)
qpt = prep_data.size > 0
# Get fitter kwargs
fitter_kwargs = {}
if meas_basis:
fitter_kwargs["measurement_basis"] = meas_basis
if meas_qubits:
fitter_kwargs["measurement_qubits"] = meas_qubits
if cond_meas_indices:
fitter_kwargs["conditional_measurement_indices"] = cond_meas_indices
if prep_basis:
fitter_kwargs["preparation_basis"] = prep_basis
if prep_qubits:
fitter_kwargs["preparation_qubits"] = prep_qubits
if cond_prep_indices:
fitter_kwargs["conditional_preparation_indices"] = cond_prep_indices
fitter_kwargs.update(**self.options.fitter_options)
fitter = self._get_fitter(self.options.fitter)
# Fit state results
state_results = self._fit_state_results(
fitter,
outcome_data,
shot_data,
meas_data,
prep_data,
qpt=qpt,
**fitter_kwargs,
)
other_results = []
# Compute fidelity with target
if len(state_results) == 1:
other_results += self._fidelity_result(
state_results[0],
fitter,
outcome_data,
shot_data,
meas_data,
prep_data,
qpt=qpt,
**fitter_kwargs,
)
# Check positive
other_results += self._positivity_result(state_results, qpt=qpt)
# Check trace preserving
if qpt:
output_dim = np.prod(state_results[0].value.output_dims())
other_results += self._tp_result(state_results, output_dim)
# Finally format state result metadata to remove eigenvectors
# which are no longer needed to reduce size
for state_result in state_results:
state_result.extra.pop("eigvecs")
analysis_results = state_results + other_results
if self.options.extra:
for res in analysis_results:
res.extra.update(self.options.extra)
return analysis_results, []
def _fit_state_results(
self,
fitter: Callable,
outcome_data: np.ndarray,
shot_data: np.ndarray,
measurement_data: np.ndarray,
preparation_data: np.ndarray,
qpt: Union[bool, str, None] = "auto",
**fitter_kwargs,
):
"""Fit state results from tomography data,"""
try:
fits, fitter_metadata = fitter(
outcome_data,
shot_data,
measurement_data,
preparation_data,
**fitter_kwargs,
)
except AnalysisError as ex:
raise AnalysisError(f"Tomography fitter failed with error: {str(ex)}") from ex
# Post process fit
states, states_metadata = postprocess_fitter(
fits,
fitter_metadata,
make_positive=self.options.rescale_positive,
trace="auto" if self.options.rescale_trace else None,
qpt=qpt,
)
# Convert to results
state_results = [
AnalysisResultData("state", state, extra=extra)
for state, extra in zip(states, states_metadata)
]
return state_results
def _fidelity_result(
self,
state_result: AnalysisResultData,
fitter: Callable,
outcome_data: np.ndarray,
shot_data: np.ndarray,
measurement_data: np.ndarray,
preparation_data: np.ndarray,
qpt: bool = False,
**fitter_kwargs,
) -> List[AnalysisResultData]:
"""Calculate fidelity result if a target has been set"""
target = self.options.target
if target is None:
return []
# Compute fidelity
name = "process_fidelity" if qpt else "state_fidelity"
fidelity = self._compute_fidelity(state_result, target, qpt=qpt)
if not self.options.target_bootstrap_samples:
# No bootstrapping
return [AnalysisResultData(name, fidelity)]
# Optionally, Estimate std error of fidelity via boostrapping
seed = self.options.target_bootstrap_seed
if isinstance(seed, np.random.Generator):
rng = seed
else:
rng = np.random.default_rng(seed)
prob_data = outcome_data / shot_data[None, :, None]
bs_fidelities = []
for _ in range(self.options.target_bootstrap_samples):
# TODO: remove conditional once numpy is pinned at 1.22 and above
if numpy_version() >= (1, 22):
sampled_data = rng.multinomial(shot_data, prob_data)
else:
sampled_data = np.zeros_like(outcome_data)
for i in range(prob_data.shape[0]):
for j in range(prob_data.shape[1]):
sampled_data[i, j] = rng.multinomial(shot_data[j], prob_data[i, j])
try:
state_results = self._fit_state_results(
fitter,
sampled_data,
shot_data,
measurement_data,
preparation_data,
qpt=qpt,
**fitter_kwargs,
)
bs_fidelities.append(self._compute_fidelity(state_results[0], target, qpt=qpt))
except AnalysisError:
pass
bs_stderr = np.std(bs_fidelities)
return [
AnalysisResultData(
name,
ufloat(fidelity, bs_stderr),
extra={"bootstrap_samples": bs_fidelities},
)
]
@staticmethod
def _positivity_result(
state_results: List[AnalysisResultData], qpt: bool = False
) -> List[AnalysisResultData]:
"""Check if eigenvalues are positive"""
total_cond = defaultdict(float)
comps_cond = defaultdict(list)
comps_pos = defaultdict(list)
name = "completely_positive" if qpt else "positive"
for result in state_results:
cond_idx = result.extra.get("conditional_measurement_index", None)
evals = result.extra["eigvals"]
# Check if component is positive and add to extra if so
cond = np.sum(np.abs(evals[evals < 0]))
pos = bool(np.isclose(cond, 0))
result.extra[name] = pos
# Add component to combined result
comps_cond[cond_idx].append(cond)
comps_pos[cond_idx].append(pos)
total_cond[cond_idx] += cond * result.extra["conditional_probability"]
# Check if combined conditional state is positive
results = []
for key, delta in total_cond.items():
is_pos = bool(np.isclose(delta, 0))
result = AnalysisResultData(name, is_pos)
if not is_pos:
result.extra = {
"delta": delta,
"components": comps_pos[key],
"components_delta": comps_cond[key],
}
if key:
result.extra["conditional_measurement_index"] = key
results.append(result)
return results
@staticmethod
def _tp_result(
state_results: List[AnalysisResultData],
output_dim: int = 1,
) -> List[AnalysisResultData]:
"""Check if QPT channel is trace preserving"""
# Construct the Kraus TP condition matrix sum_i K_i^dag K_i
# summed over all components k
kraus_cond = {}
for result in state_results:
evals = result.extra["eigvals"]
evecs = result.extra["eigvecs"]
prob = result.extra["conditional_probability"]
cond_meas_idx = result.extra.get("conditional_measurement_index", None)
cond_prep_idx = result.extra.get("conditional_preparation_index", None)
cond_idx = (cond_prep_idx, cond_meas_idx)
size = len(evals)
input_dim = size // output_dim
mats = np.reshape(evecs.T, (size, input_dim, output_dim), order="F")
comp_cond = np.einsum("i,iaj,ibj->ab", evals, mats.conj(), mats)
if cond_idx in kraus_cond:
kraus_cond[cond_idx] += prob * comp_cond
else:
kraus_cond[cond_idx] = prob * comp_cond
results = []
for key, val in kraus_cond.items():
tp_cond = np.sum(np.abs(la.eigvalsh(val - np.eye(input_dim))))
is_tp = bool(np.isclose(tp_cond, 0, atol=1e-5))
result = AnalysisResultData("trace_preserving", is_tp, extra={})
if not is_tp:
result.extra["delta"] = tp_cond
if key:
result.extra["conditional_measurement_index"] = key
results.append(result)
return results
@staticmethod
def _compute_fidelity(
state_result: AnalysisResultData,
target: Union[Choi, DensityMatrix],
qpt: bool = False,
) -> AnalysisResultData:
"""Faster computation of fidelity from eigen decomposition"""
if qpt:
input_dim = np.prod(state_result.value.input_dims())
else:
input_dim = 1
evals = state_result.extra["eigvals"]
evecs = state_result.extra["eigvecs"]
# Format target to statevector or densitymatrix array
if target is None:
raise AnalysisError("No target state provided")
if isinstance(target, QuantumChannel):
target_state = Choi(target).data / input_dim
elif isinstance(target, BaseOperator):
target_state = np.ravel(Operator(target), order="F") / np.sqrt(input_dim)
else:
# Statevector or density matrix
target_state = np.array(target)
if target_state.ndim == 1:
rho = evecs @ (evals / input_dim * evecs).T.conj()
fidelity = np.real(target_state.conj() @ rho @ target_state)
else:
sqrt_rho = evecs @ (np.sqrt(evals / input_dim) * evecs).T.conj()
eig = la.eigvalsh(sqrt_rho @ target_state @ sqrt_rho)
fidelity = np.sum(np.sqrt(np.maximum(eig, 0))) ** 2
return fidelity