Skip to content

Commit

Permalink
Update tomography analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
@chriseclectic
chriseclectic committed Feb 13, 2023
1 parent e1c0295 commit d360b8c
Showing 1 changed file with 34 additions and 262 deletions.
296 changes: 34 additions & 262 deletions qiskit_experiments/library/tomography/tomography_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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,
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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

0 comments on commit d360b8c

Please sign in to comment.