From d360b8c3d3eea85c989d8d143c572bafc5cb2af5 Mon Sep 17 00:00:00 2001 From: "Christopher J. Wood" Date: Mon, 13 Feb 2023 18:09:52 -0500 Subject: [PATCH] Update tomography analysis --- .../library/tomography/tomography_analysis.py | 296 ++---------------- 1 file changed, 34 insertions(+), 262 deletions(-) diff --git a/qiskit_experiments/library/tomography/tomography_analysis.py b/qiskit_experiments/library/tomography/tomography_analysis.py index e1529a0d04..678e468dc9 100644 --- a/qiskit_experiments/library/tomography/tomography_analysis.py +++ b/qiskit_experiments/library/tomography/tomography_analysis.py @@ -14,9 +14,8 @@ """ -from typing import List, Dict, Tuple, Union, Optional, Callable +from typing import List, Union, Callable import warnings -import time import numpy as np import scipy.linalg as la @@ -26,9 +25,9 @@ from qiskit_experiments.exceptions import AnalysisError from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options -from .basis import MeasurementBasis, PreparationBasis from .fitters import ( tomography_fitter_data, + postprocess_fitter, linear_inversion, scipy_linear_lstsq, scipy_gaussian_lstsq, @@ -153,254 +152,73 @@ def _run_analysis(self, experiment_data): if measurement_basis and measurement_qubits: outcome_shape = measurement_basis.outcome_shape(measurement_qubits) - fitter_data = tomography_fitter_data( + outcome_data, shot_data, meas_data, prep_data = tomography_fitter_data( experiment_data.data(), outcome_shape=outcome_shape, ) + fitter = self._get_fitter(self.options.fitter) - # Run fitter - fits, fitter_metadata = self._run_fit( - self._get_fitter(self.options.fitter), - *fitter_data, - measurement_basis=measurement_basis, - measurement_qubits=measurement_qubits, - preparation_basis=preparation_basis, - preparation_qubits=preparation_qubits, - conditional_measurement_indices=conditional_measurement_indices, - fitter_options=self.options.fitter_options, - ) - - # Post process fit - analysis_results = self._postprocess_fits( - fits, - fitter_metadata=fitter_metadata, - make_positive=self.options.rescale_positive, - target_state=self.options.target, - ) - return analysis_results, [] - - def _run_fit( - self, - fitter: Callable, - outcome_data: np.ndarray, - shot_data: np.ndarray, - measurement_data: np.ndarray, - preparation_data: np.ndarray, - measurement_basis: Optional[MeasurementBasis] = None, - preparation_basis: Optional[PreparationBasis] = None, - measurement_qubits: Optional[List[int]] = None, - preparation_qubits: Optional[List[int]] = None, - conditional_measurement_indices: Optional[List[int]] = None, - fitter_options: dict = None, - ): - """Run tomography fitter on tomography data""" - # Construct default values for qubit options if not provided - if preparation_qubits is None: - preparation_qubits = tuple(range(preparation_data.shape[1])) - if measurement_qubits is None: - measurement_qubits = tuple(range(measurement_data.shape[1])) - - # Get dimension of the preparation and measurement qubits subsystems - prep_dims = (1,) - if preparation_qubits: - if not preparation_basis: - raise AnalysisError("No tomography preparation basis provided.") - prep_dims = preparation_basis.matrix_shape(preparation_qubits) - meas_dims = (1,) - full_meas_qubits = measurement_qubits - if measurement_qubits: - if conditional_measurement_indices is not None: - # Remove conditional qubits from full meas qubits - full_meas_qubits = [ - q - for i, q in enumerate(measurement_qubits) - if i not in conditional_measurement_indices - ] - if full_meas_qubits: - if not measurement_basis: - raise AnalysisError("No tomography measurement basis provided.") - meas_dims = measurement_basis.matrix_shape(full_meas_qubits) - - if full_meas_qubits: - # QPT or QST - input_dims = prep_dims - output_dims = meas_dims - else: - # QST of POVM effects - input_dims = meas_dims - output_dims = prep_dims - - # Use preparation dim to set the expected trace of the fitted state. - # For QPT this is the input dimension, for QST this will always be 1. - trace = np.prod(prep_dims) if self.options.rescale_trace else None - - # Get tomography fitter options - if fitter_options is None: - fitter_options = {} - - # Work around to set proper trace and trace preserving constraints for - # cvxpy fitter - if fitter in self._cvxpy_fitters: - fitter_options = fitter_options.copy() - - # Add default value for CVXPY trace constraint if no user value is provided - # Use preparation dim to set the expected trace of the fitted state. - # For QPT this is the input dimension, for QST this will always be 1. - if "trace" not in fitter_options: - fitter_options["trace"] = trace - - # By default add trace preserving constraint to cvxpy QPT fit - if preparation_data.shape[1] > 0 and "trace_preserving" not in fitter_options: - fitter_options["trace_preserving"] = True - - # Run tomography fitter - t_fitter_start = time.time() try: - fit, fitter_metadata = fitter( + fits, fitter_metadata = fitter( outcome_data, shot_data, - measurement_data, - preparation_data, + meas_data, + prep_data, measurement_basis=measurement_basis, - preparation_basis=preparation_basis, measurement_qubits=measurement_qubits, + preparation_basis=preparation_basis, preparation_qubits=preparation_qubits, conditional_measurement_indices=conditional_measurement_indices, - **fitter_options, + **self.options.fitter_options, ) except AnalysisError as ex: raise AnalysisError(f"Tomography fitter failed with error: {str(ex)}") from ex - t_fitter_stop = time.time() - - # Add fitter metadata - if fitter_metadata is None: - fitter_metadata = {} - fitter_metadata["fitter"] = fitter.__name__ - fitter_metadata["fitter_time"] = t_fitter_stop - t_fitter_start - fitter_metadata["input_dims"] = input_dims - fitter_metadata["output_dims"] = output_dims - if trace is not None: - fitter_metadata["trace"] = trace - return fit, fitter_metadata - @classmethod - def _postprocess_fits( - cls, - fits: List[np.ndarray], - fitter_metadata: Optional[Dict] = None, - make_positive: bool = False, - target_state: Optional[Union[Choi, DensityMatrix]] = None, - ) -> Dict[str, any]: - """Post-process raw fitter result""" - # Get dimension and trace from fitter metadata - trace = fitter_metadata.get("trace", None) - conditionals = fitter_metadata.pop("component_conditionals", None) - input_dims = fitter_metadata.get("input_dims", None) - output_dims = fitter_metadata.get("output_dims", None) - - # Convert fitter matrix to state data for post-processing - input_dim = np.prod(input_dims) if input_dims else 1 - qpt = input_dim > 1 - state_results = [ - cls._state_result( - fit, - make_positive=make_positive, - trace=trace, - input_dims=input_dims, - output_dims=output_dims, - fitter_metadata=fitter_metadata, - ) - for fit in fits - ] + # 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, + ) - # Compute the conditional probability of each component so that the - # total probability of all components is 1, and optional rescale trace - # of each component - fit_traces = [res.extra.pop("fit_trace") for res in state_results] - total_trace = sum(fit_traces) - for i, (fit_trace, res) in enumerate(zip(fit_traces, state_results)): - # Compute conditional component probability from the the component - # non-rescaled fit trace - res.extra["component_probability"] = fit_trace / total_trace - res.extra["component_index"] = i - if conditionals: - res.extra["component_conditional"] = conditionals[i] + if prep_data.shape[1]: + qpt = True + input_dim = np.prod(states[0].input_dims()) + else: + qpt = False + input_dim = 1 + # Convert to results + state_results = [ + AnalysisResultData("state", state, extra=extra) + for state, extra in zip(states, states_metadata) + ] other_results = [] + # Compute fidelity with target + target_state = self.options.target if len(state_results) == 1 and target_state is not None: # Note: this currently only works for non-conditional tomography other_results.append( - cls._fidelity_result(state_results[0], target_state, input_dim=input_dim) + self._fidelity_result(state_results[0], target_state, input_dim=input_dim) ) # Check positive - other_results.append(cls._positivity_result(state_results, qpt=qpt)) + other_results.append(self._positivity_result(state_results, qpt=qpt)) # Check trace preserving if qpt: - other_results.append(cls._tp_result(state_results, input_dim=input_dim)) + other_results.append(self._tp_result(state_results, input_dim=input_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") - return state_results + other_results - - @classmethod - def _state_result( - cls, - fit: np.ndarray, - make_positive: bool = False, - trace: Optional[float] = None, - input_dims: Optional[Tuple[int, ...]] = None, - output_dims: Optional[Tuple[int, ...]] = None, - fitter_metadata: Optional[Dict] = None, - ) -> List[AnalysisResultData]: - """Convert fit data to state result data""" - # Get eigensystem of state fit - raw_eigvals, eigvecs = cls._state_eigensystem(fit) - - # Optionally rescale eigenvalues to be non-negative - if make_positive and np.any(raw_eigvals < 0): - eigvals = cls._make_positive(raw_eigvals) - fit = eigvecs @ (eigvals * eigvecs).T.conj() - rescaled_psd = True - else: - eigvals = raw_eigvals - rescaled_psd = False - - # Optionally rescale fit trace - fit_trace = np.sum(eigvals) - if ( - trace is not None - and not np.isclose(fit_trace, 0, atol=1e-10) - and not np.isclose(abs(fit_trace - trace), 0, atol=1e-10) - ): - scale = trace / fit_trace - fit = fit * scale - eigvals = eigvals * scale - else: - trace = fit_trace + analysis_results = state_results + other_results - # Convert class of value - if input_dims and np.prod(input_dims) > 1: - value = Choi(fit, input_dims=input_dims, output_dims=output_dims) - else: - value = DensityMatrix(fit, dims=output_dims) - - # Construct state result extra metadata - extra = { - "trace": trace, - "eigvals": eigvals, - "raw_eigvals": raw_eigvals, - "rescaled_psd": rescaled_psd, - "fit_trace": fit_trace, - "eigvecs": eigvecs, - "fitter_metadata": fitter_metadata or {}, - } - return AnalysisResultData("state", value, extra=extra) + return analysis_results, [] @staticmethod def _positivity_result( @@ -491,49 +309,3 @@ def _fidelity_result( eig = la.eigvalsh(sqrt_rho @ target_state @ sqrt_rho) fidelity = np.sum(np.sqrt(np.maximum(eig, 0))) ** 2 return AnalysisResultData(name, fidelity) - - @staticmethod - def _state_eigensystem(fit: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: - """Compute the eigensystem of the fitted state. - - The eigenvalues are returned as a real array ordered from - smallest to largest eigenvalues. - - Args: - fit: the fitted state matrix. - - Returns: - A pair of (eigenvalues, eigenvectors). - """ - evals, evecs = la.eigh(fit) - # Truncate eigenvalues to real part - evals = np.real(evals) - # Sort eigensystem from largest to smallest eigenvalues - sort_inds = np.flip(np.argsort(evals)) - return evals[sort_inds], evecs[:, sort_inds] - - @staticmethod - def _make_positive(evals: np.ndarray, epsilon: float = 0) -> np.ndarray: - """Rescale a real vector to be non-negative. - - This truncates any negative values to zero and rescales - the remaining eigenvectors such that the sum of the vector - is preserved. - """ - if epsilon < 0: - raise AnalysisError("epsilon must be non-negative.") - ret = evals.copy() - dim = len(evals) - idx = dim - 1 - accum = 0.0 - while idx >= 0: - shift = accum / (idx + 1) - if evals[idx] + shift < epsilon: - ret[idx] = 0 - accum = accum + evals[idx] - idx -= 1 - else: - for j in range(idx + 1): - ret[j] = evals[j] + shift - break - return ret