Bugfix reparameterize

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 <ideal|exp(errorgen) = <effect|")
+                    errgen_vec = _np.linalg.pinv(phys_directions)  @ soln.x
+                    errorgen.from_vector(errgen_vec)
+
                 EffectiveExpErrorgen = _IdentityPlusErrorgenOp if lndtype.meta == '1+' else _ExpErrorgenOp
                 return ComposedPOVM(EffectiveExpErrorgen(errorgen), base_povm, mx_basis=basis)
 

pygsti/modelmembers/states/__init__.py

@@ -32,6 +32,7 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.tools import basistools as _bt
 from pygsti.tools import optools as _ot
+from pygsti.tools import sum_of_negative_choi_eigenvalues_gate
 
 # Avoid circular import
 import pygsti.modelmembers as _mm
@@ -258,20 +259,54 @@ def convert(state, to_type, basis, ideal_state=None, flatten_structure=False):
                 errorgen = _LindbladErrorgen.from_error_generator(state.state_space.dim, to_type, proj_basis,
                                                                   basis, truncate=True, evotype=state.evotype)
                 if st is not state and not _np.allclose(st.to_dense(), state.to_dense()):
-                    #Need to set errorgen so exp(errorgen)|st> == |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))

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)
 

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):
     """

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