Bugfix reparameterize

pygsti/modelmembers/povms/__init__.py

@@ -11,8 +11,10 @@
 #***************************************************************************************************
 import functools as _functools
 import itertools as _itertools
-
+import warnings
 import numpy as _np
+import scipy.linalg as _spl
+import scipy.optimize as _spo
 
 from .complementeffect import ComplementPOVMEffect
 from .composedeffect import ComposedPOVMEffect
@@ -34,6 +36,8 @@
 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
+from pygsti.baseobjs import Basis
 
 # Avoid circular import
 import pygsti.modelmembers as _mm
@@ -418,7 +422,7 @@ def povm_type_from_op_type(op_type):
     return povm_type_preferences
 
 
-def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False):
+def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False, cp_penalty=1e-7):
     """
     TODO: update docstring
     Convert a POVM to a new type of parameterization.
@@ -450,13 +454,22 @@ def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False):
         are separately converted, leaving the original POVM's structure
         unchanged.  When `True`, composed and embedded operations are "flattened"
         into a single POVM of the requested `to_type`.
+
+    cp_penalty : float, optional (default 1e-7)
+            Converting SPAM operations to an error generator representation may 
+            introduce trivial gauge degrees of freedom. These gauge degrees of freedom 
+            are called trivial because they quite literally do not change the dense representation 
+            (i.e. Hilbert-Schmidt vectors) at all. Despite being trivial, error generators along 
+            this trivial gauge orbit may be non-CP, so this cptp penalty is used to favor channels 
+            within this gauge orbit which are CPTP.
 
     Returns
     -------
     POVM
         The converted POVM vector, usually a distinct
         object from the object passed as input.
     """
+
     to_types = to_type if isinstance(to_type, (tuple, list)) else (to_type,)  # HACK to support multiple to_type values
     error_msgs = {}
 
@@ -487,6 +500,7 @@ def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False):
                 from ..operations import LindbladParameterization as _LindbladParameterization
                 lndtype = _LindbladParameterization.cast(to_type)
 
+
                 #Construct a static "base" POVM
                 if isinstance(povm, ComputationalBasisPOVM):  # special easy case
                     base_povm = ComputationalBasisPOVM(povm.state_space.num_qubits, povm.evotype)  # just copy it?
@@ -498,14 +512,105 @@ 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)
+
+                #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_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 trivial gauge", and this function is meant to avoid it
+                def calc_physical_subspace(dense_ideal_povm, epsilon = 1e-4):
+
+                    degrees_of_freedom = (dense_ideal_povm.shape[0] - 1) * dense_ideal_povm.shape[1]
+                    errgen = _LindbladErrorgen.from_error_generator(povm.state_space.dim, parameterization=to_type)
+
+                    if degrees_of_freedom > errgen.num_params:
+                        warnings.warn("POVM has more degrees of freedom than the available number of parameters, representation in this parameterization is not guaranteed")
+                    exp_errgen = _ExpErrorgenOp(errgen)
+
+                    num_errgens = errgen.num_params
+                    #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.size
+
+                    #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))
+                    new_vec = _np.zeros(num_errgens)
+                    for i in range(num_errgens):
+
+                        new_vec[i] = epsilon
+                        exp_errgen.from_vector(new_vec)
+                        new_vec[i] = 0
+                        vectorized_povm = _np.zeros(num_entries)
+                        perturbed_povm = (dense_ideal_povm @ exp_errgen.to_dense() - dense_ideal_povm)/epsilon 
+
+                        vectorized_povm = perturbed_povm.flatten(order='F')
+
+                        J[:,i] = vectorized_povm
+
+                    _,S,Vt = _np.linalg.svd(J, full_matrices=False)
+
+                    #Only return nontrivial singular vectors
+                    Vt = Vt[S > 1e-13, :].reshape((-1, Vt.shape[1]))
+                    return Vt
+
+
+                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):
+
+                    #For some reason adding the sum_of_negative_choi_eigenvalues_gate term
+                    #resulted in minimize() sometimes choosing NaN values for v. There are
+                    #two stack exchange issues showing this problem with no solution.
+                    if _np.isnan(v).any():
+                        v = _np.zeros(len(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)
+                    proc_matrix = _spl.expm(errorgen.to_dense())
+
+                    return _np.linalg.norm(dense_povm - dense_ideal_povm @ proc_matrix) + cp_penalty * sum_of_negative_choi_eigenvalues_gate(proc_matrix, basis)
+
+                soln = _spo.minimize(_objfn, _np.zeros(len(phys_directions), 'd'), method="Nelder-Mead", options={},
+                                        tol=1e-13) 
+                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.lstsq(phys_directions, soln.x)[0]
+                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

@@ -29,7 +29,9 @@
 from .tpstate import TPState
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.tools import basistools as _bt
+from pygsti.baseobjs import Basis
 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
@@ -185,7 +187,7 @@ def state_type_from_op_type(op_type):
     return state_type_preferences
 
 
-def convert(state, to_type, basis, ideal_state=None, flatten_structure=False):
+def convert(state, to_type, basis, ideal_state=None, flatten_structure=False, cp_penalty=1e-7):
     """
     TODO: update docstring
     Convert SPAM vector to a new type of parameterization.
@@ -217,6 +219,10 @@ def convert(state, to_type, basis, ideal_state=None, flatten_structure=False):
         are separately converted, leaving the original state's structure
         unchanged.  When `True`, composed and embedded operations are "flattened"
         into a single state of the requested `to_type`.
+
+    cp_penalty : float, optional
+        CPTP penalty that gets factored into the optimization to find the resulting model
+        when converting to an error-generator type.
 
     Returns
     -------
@@ -234,7 +240,6 @@ def convert(state, to_type, basis, ideal_state=None, flatten_structure=False):
                          'static unitary': StaticPureState,
                          'static clifford': ComputationalBasisState}
     NoneType = type(None)
-
     for to_type in to_types:
         try:
             if isinstance(state, destination_types.get(to_type, NoneType)):
@@ -256,20 +261,57 @@ 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
+
+
+
+                    #GLND for states suffers from "trivial gauge" freedom. This function identifies
+                    #the physical directions to avoid this gauge.
+                    def calc_physical_subspace(ideal_prep, epsilon = 1e-4):
+                        errgen = _LindbladErrorgen.from_error_generator(2**(2*num_qubits), parameterization=to_type)
+                        num_errgens = errgen.num_params
+
+                        exp_errgen = _ExpErrorgenOp(errgen)
+
+                        #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(num_errgens):
+                            new_vec = _np.zeros(num_errgens)
+                            new_vec[i] = epsilon
+                            exp_errgen.from_vector(new_vec)
+                            J[:,i] = (exp_errgen.to_dense() @ ideal_prep - ideal_prep)[1:]/epsilon
+
+                        _,S,Vt = _np.linalg.svd(J, full_matrices=False)
+
+                        #Only return nontrivial singular vectors
+                        Vt = Vt[S > 1e-13, :].reshape((-1, Vt.shape[1]))
+                        return Vt
+
+                    phys_directions = calc_physical_subspace(dense_state)
+
+                    #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)
-                    #print("DEBUG: opt done: ",soln.success, soln.fun, soln.x)  # REMOVE
+                        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) + cp_penalty * sum_of_negative_choi_eigenvalues_gate(proc_matrix, basis)
+
+                    soln = _spo.minimize(_objfn, _np.zeros(len(phys_directions), 'd'), method="Nelder-Mead", options={},
+                                         tol=1e-13)  
+
                     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.lstsq(phys_directions, soln.x)[0]
+                    errorgen.from_vector(errgen_vec)
 
                 EffectiveExpErrorgen = _IdentityPlusErrorgenOp if lndtype.meta == '1+' else _ExpErrorgenOp
                 return ComposedState(st, EffectiveExpErrorgen(errorgen))
@@ -283,7 +325,6 @@ def _objfn(v):
                     return create_from_dmvec(vec, to_type, basis, state.evotype, state.state_space)
 
         except ValueError as e:
-            #_warnings.warn('Failed to convert state to type %s with error: %s' % (to_type, e))
             error_msgs[to_type] = str(e)  # try next to_type
 
     raise ValueError("Could not convert state to to type(s): %s\n%s" % (str(to_types), str(error_msgs)))