diff --git a/pygsti/modelmembers/povms/__init__.py b/pygsti/modelmembers/povms/__init__.py index e946abd51..b020aef84 100644 --- a/pygsti/modelmembers/povms/__init__.py +++ b/pygsti/modelmembers/povms/__init__.py @@ -14,6 +14,8 @@ import itertools as _itertools import numpy as _np +import scipy.linalg as _spl +import scipy.optimize as _spo from .complementeffect import ComplementPOVMEffect from .composedeffect import ComposedPOVMEffect @@ -334,6 +336,8 @@ def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False): The converted POVM vector, usually a distinct object from the object passed as input. """ + + ##TEST CONVERSION BETWEEN LINBLAD TYPES to_types = to_type if isinstance(to_type, (tuple, list)) else (to_type,) # HACK to support multiple to_type values error_msgs = {} @@ -375,14 +379,96 @@ def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False): for lbl, vec in povm.items()] else: raise RuntimeError('Evotype must be compatible with Hilbert ops to use pure effects') - except Exception: # try static mixed states next: - base_items = [(lbl, convert_effect(vec, 'static', basis, idl.get(lbl, None), flatten_structure)) - for lbl, vec in povm.items()] + except RuntimeError: # try static mixed states next: + #if idl.get(lbl,None) is not None: + + base_items = [] + for lbl, vec in povm.items(): + ideal_effect = idl.get(lbl,None) + if ideal_effect is not None: + base_items.append((lbl, convert_effect(ideal_effect, 'static', basis, ideal_effect, flatten_structure))) + else: + base_items.append((lbl, convert_effect(vec, 'static', basis, idl.get(lbl, None), flatten_structure))) base_povm = UnconstrainedPOVM(base_items, povm.evotype, povm.state_space) proj_basis = 'PP' if povm.state_space.is_entirely_qubits else basis errorgen = _LindbladErrorgen.from_error_generator(povm.state_space.dim, lndtype, proj_basis, basis, truncate=True, evotype=povm.evotype) + if to_type == 'GLND' and isinstance(povm, destination_types.get('full TP', NoneType)): + + #Collect all ideal effects + base_dense_effects = [] + for item in base_items: + dense_effect = item[1].to_dense() + base_dense_effects.append(dense_effect.reshape((1,len(dense_effect)))) + + dense_ideal_povm = _np.concatenate(base_dense_effects, axis=0) + + #Collect all noisy effects + dense_effects = [] + for effect in povm.values(): + dense_effect = effect.to_dense() + dense_effects.append(dense_effect.reshape((1,len(dense_effect)))) + + dense_ideal_povm = _np.concatenate(base_dense_effects, axis=0) + dense_povm = _np.concatenate(dense_effects, axis=0) + + #It is often the case that there are more error generators than physical degrees of freedom in the POVM + #We define a function which finds linear comb. of errgens that span these degrees of freedom. + #This has been called "the dumb gauge", and this function is meant to avoid it + def calc_physical_subspace(dense_ideal_povm, epsilon = 1e-9): + + errgen = _LindbladErrorgen.from_error_generator(4, parameterization="GLND") + exp_errgen = _ExpErrorgenOp(errgen) + num_qubits = povm.state_space.num_qubits + num_errgens = 2**(4*num_qubits)-2**(2*num_qubits) + #TODO: Maybe we can use the num of params instead of number of matrix entries, as some of them are linearly dependent. + #i.e E0 completely determines E1 if those are the only two povm elements (E0 + E1 = Identity) + num_entries = dense_ideal_povm.shape[0]*dense_ideal_povm.shape[1] + assert num_errgens > povm.num_params, "POVM has too many elements, GLND parameterization is not possible" + + ideal_vec = _np.zeros(num_errgens) + + #Compute the jacobian with respect to the error generators. This will allow us to see which + #error generators change the POVM entries + J = _np.zeros((num_entries,num_errgens)) + + for i in range(len(ideal_vec)): + new_vec = ideal_vec.copy() + new_vec[i] = epsilon + exp_errgen.from_vector(new_vec) + vectorized_povm = _np.zeros(num_entries) + perturbed_povm = (dense_ideal_povm @ exp_errgen.to_dense() - dense_ideal_povm)/epsilon + + perturbed_povm_t = perturbed_povm.transpose() + for j, column in enumerate(perturbed_povm_t): + vectorized_povm[j*len(perturbed_povm_t[0]):(j+1)*len(perturbed_povm_t[0])] = column + + J[:,i] = vectorized_povm.transpose() + + _,S,V = _np.linalg.svd(J) + return V[:len(S),] + + phys_directions = calc_physical_subspace(dense_ideal_povm) + + #We use optimization to find the best error generator representation + #we only vary physical directions, not independent error generators + def _objfn(v): + L_vec = _np.zeros(len(phys_directions[0])) + for coeff, phys_direction in zip(v,phys_directions): + L_vec += coeff * phys_direction + errorgen.from_vector(L_vec) + return _np.linalg.norm(dense_povm - dense_ideal_povm @ _spl.expm(errorgen.to_dense())) + + #def callback(x): print("callbk: ",_np.linalg.norm(x),_objfn(x)) # REMOVE + soln = _spo.minimize(_objfn, _np.zeros(len(phys_directions), 'd'), method="Nelder-Mead", options={}, + tol=1e-13) # , callback=callback) + #print("DEBUG: opt done: ",soln.success, soln.fun, soln.x) # REMOVE + if not soln.success and soln.fun > 1e-6: # not "or" because success is often not set correctly + raise ValueError("Failed to find an errorgen such that == |state> + dense_st = st.to_dense() dense_state = state.to_dense() + num_qubits = st.state_space.num_qubits + num_errgens = 2**(4*num_qubits)-2**(2*num_qubits) + + #GLND for states suffers from "dumb gauge" freedom. This function identifies + #the physical directions to avoid this gauge. + def calc_physical_subspace(ideal_prep, epsilon = 1e-9): + + errgen = _LindbladErrorgen.from_error_generator(2**(2*num_qubits), parameterization="GLND") + exp_errgen = _ExpErrorgenOp(errgen) + ideal_vec = _np.zeros(num_errgens) + + #Compute the jacobian with respect to the error generators. This will allow us to see which + #error generators change the POVM entries + J = _np.zeros((state.num_params, num_errgens)) + + for i in range(len(ideal_vec)): + new_vec = ideal_vec.copy() + new_vec[i] = epsilon + exp_errgen.from_vector(new_vec) + J[:,i] = (exp_errgen.to_dense() @ ideal_prep - ideal_prep)[1:]/epsilon + + _,S,V = _np.linalg.svd(J) + return V[:len(S),] + + phys_directions = calc_physical_subspace(dense_state) + cptp_penalty = .5 + #We use optimization to find the best error generator representation + #we only vary physical directions, not independent error generators def _objfn(v): - errorgen.from_vector(v) - return _np.linalg.norm(_spl.expm(errorgen.to_dense()) @ dense_st - dense_state) - #def callback(x): print("callbk: ",_np.linalg.norm(x),_objfn(x)) # REMOVE - soln = _spo.minimize(_objfn, _np.zeros(errorgen.num_params, 'd'), method="CG", options={}, - tol=1e-8) # , callback=callback) + L_vec = _np.zeros(len(phys_directions[0])) + for coeff, phys_direction in zip(v,phys_directions): + L_vec += coeff * phys_direction + errorgen.from_vector(L_vec) + proc_matrix = _spl.expm(errorgen.to_dense()) + return _np.linalg.norm(proc_matrix @ dense_st - dense_state) + cptp_pentaly * sum_of_negative_choi_eigenvalues_gate(proc_matrix) + + soln = _spo.minimize(_objfn, _np.zeros(len(phys_directions), 'd'), method="Nelder-Mead", options={}, + tol=1e-13) # , callback=callback) #print("DEBUG: opt done: ",soln.success, soln.fun, soln.x) # REMOVE if not soln.success and soln.fun > 1e-6: # not "or" because success is often not set correctly raise ValueError("Failed to find an errorgen such that exp(errorgen)|ideal> = |state>") - errorgen.from_vector(soln.x) + + errgen_vec = _np.linalg.pinv(phys_directions) @ soln.x + errorgen.from_vector(errgen_vec) EffectiveExpErrorgen = _IdentityPlusErrorgenOp if lndtype.meta == '1+' else _ExpErrorgenOp return ComposedState(st, EffectiveExpErrorgen(errorgen)) diff --git a/pygsti/models/modelparaminterposer.py b/pygsti/models/modelparaminterposer.py index f92f56a70..3d0ea4f75 100644 --- a/pygsti/models/modelparaminterposer.py +++ b/pygsti/models/modelparaminterposer.py @@ -66,7 +66,7 @@ def __init__(self, transform_matrix): self.transform_matrix = transform_matrix # cols specify a model parameter in terms of op params. self.inv_transform_matrix = _np.linalg.pinv(self.transform_matrix) super().__init__(transform_matrix.shape[1], transform_matrix.shape[0]) - + def model_paramvec_to_ops_paramvec(self, v): return _np.dot(self.transform_matrix, v) diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py index e204c5206..055b65543 100644 --- a/pygsti/tools/jamiolkowski.py +++ b/pygsti/tools/jamiolkowski.py @@ -273,6 +273,27 @@ def fast_jamiolkowski_iso_std_inv(choi_mx, op_mx_basis): #transform operation matrix into appropriate basis return _bt.change_basis(opMxInStdBasis, op_mx_basis.create_equivalent('std'), op_mx_basis) +def sum_of_negative_choi_eigenvalues_gate(op_mx, op_mx_basis): + """ + Compute the sum of the negative Choi eigenvalues of a process matrix. + + Parameters + ---------- + op_mx : np.array + + op_mx_basis : Basis + + Returns + ------- + float + the sum of the negative eigenvalues of the Choi representation of op_mx + """ + sumOfNeg = 0 + J = fast_jamiolkowski_iso_std(op_mx, op_mx_basis) # Choi mx basis doesn't matter + evals = _np.linalg.eigvals(J) # could use eigvalsh, but wary of this since eigh can be wrong... + for ev in evals: + if ev.real < 0: sumOfNeg -= ev.real + return sumOfNeg def sum_of_negative_choi_eigenvalues(model, weights=None): """ diff --git a/test/unit/objects/test_model.py b/test/unit/objects/test_model.py index fdad0a22c..cf5f24460 100644 --- a/test/unit/objects/test_model.py +++ b/test/unit/objects/test_model.py @@ -136,6 +136,16 @@ def test_set_all_parameterizations_HS(self): nParamsPerGate=6, nParamsPerSP=6 ) + + def test_set_all_parameterizations_(self): + self.model.set_all_parameterizations("H+S") + self._assert_model_params( + nOperations=3, + nSPVecs=2, + nEVecs=0, + nParamsPerGate=6, + nParamsPerSP=6 + ) def test_element_accessors(self): # XXX what does this test cover and is it useful? EGN: covers the __getitem__/__setitem__ functions of model