WIP: CPTP instruments

pygsti/baseobjs/__init__.py

@@ -13,7 +13,7 @@
 from .smartcache import SmartCache
 from .verbosityprinter import VerbosityPrinter
 from .profiler import Profiler
-from .basis import Basis, BuiltinBasis, ExplicitBasis, TensorProdBasis, DirectSumBasis
+from .basis import Basis, BuiltinBasis, ExplicitBasis, TensorProdBasis, DirectSumBasis, BasisLike
 from .label import Label, CircuitLabel
 from .nicelyserializable import NicelySerializable
 from .mongoserializable import MongoSerializable

pygsti/modelmembers/instruments/__init__.py

@@ -13,8 +13,18 @@
 from .instrument import Instrument
 from .tpinstrument import TPInstrument
 from .tpinstrumentop import TPInstrumentOp
+from ..operations.cptrop import RootConjOperator, SummedOperator
 
-from pygsti.tools import optools as _ot
+import warnings as _warnings
+import scipy.linalg as _la
+import numpy as _np
+from pygsti.tools import optools as _ot, basistools as _bt
+from types import NoneType
+from pygsti.baseobjs.label import Label
+from pygsti.baseobjs.basis import Basis
+from pygsti.modelmembers import operations as _op
+from pygsti.modelmembers import povms as _pv
+from pygsti.modelmembers.povms.basepovm import _BasePOVM
 
 # Avoid circular import
 import pygsti.modelmembers as _mm
@@ -37,27 +47,37 @@ def instrument_type_from_op_type(op_type):
 
     # Limited set (only matching what is in convert)
     instr_conversion = {
-        'auto': 'full',
+        'auto': 'full TP',
         'static unitary': 'static unitary',
         'static clifford': 'static clifford',
         'static': 'static',
         'full': 'full',
         'full TP': 'full TP',
         'full CPTP': 'full CPTP',
         'full unitary': 'full unitary',
+        'GLND': 'full TP',
+        # ^ It's pretty harmless to associate GLND operations with "full TP"
+        #   instruments. In both cases we're relaxing positivity constraints.
+        'CPTPLND': 'CPTPLND'
     }
 
     instr_type_preferences = []
     for typ in op_type_preferences:
-        instr_type = None
-        if _ot.is_valid_lindblad_paramtype(typ):
-            # Lindblad types are passed through as TP only (matching current convert logic)
-            instr_type = "full TP"
-        else:
-            instr_type = instr_conversion.get(typ, None)
+        instr_type = instr_conversion.get(typ, None)
 
-        if instr_type is None:
-            continue
+        if instr_type is None and _ot.is_valid_lindblad_paramtype(typ):
+            # NOTE: need to update the message below if more lindblad
+            # types are added as keys to the instr_conversion dict.
+            msg = \
+            f"""
+            Operation type {typ} is a Lindblad parameterization, but
+            is neither 'GLND' or 'CPTPLND'. That means you might be
+            asking for a reduced-order model that's a subset of all
+            CPTP models. We don't support that parameterization at this
+            time. We're falling back to a 'full TP' parameterization!
+            """
+            _warnings.warn(msg)
+            instr_type = 'full TP' # non-CPTPLND falls back to full TP.
 
         if instr_type not in instr_type_preferences:
             instr_type_preferences.append(instr_type)
@@ -110,26 +130,153 @@ def convert(instrument, to_type, basis, ideal_instrument=None, flatten_structure
         The converted instrument, 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
+    if not isinstance(to_type, str):
+        if len(to_type) > 1:
+            raise ValueError(f"Expected to_type to be a string, but got {to_type}")
+        to_type = to_type[0]
+        assert isinstance(to_type, str)
+
     destination_types = {'full TP': TPInstrument}
-    NoneType = type(None)
-
-    for to_type in to_types:
-        try:
-            if isinstance(instrument, destination_types.get(to_type, NoneType)):
-                return instrument
-
-            if to_type == "full TP":
-                return TPInstrument(list(instrument.items()), instrument.evotype, instrument.state_space)
-            elif to_type in ("full", "static", "static unitary"):
-                from ..operations import convert as _op_convert
-                ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {}
-                members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure))
-                           for k, g in instrument.items()]
-                return Instrument(members, instrument.evotype, instrument.state_space)
+
+    if isinstance(instrument, destination_types.get( to_type, NoneType ) ):
+        return instrument
+
+    if to_type == "full TP":
+        inst_arrays = dict()
+        for k, v in instrument.items():
+            if hasattr(v, 'to_dense'):
+                inst_arrays[k] = v.to_dense('HilbertSchmidt')
             else:
-                raise ValueError("Cannot convert an instrument to type %s" % to_type)
-        except:
-            pass  # try next to_type
+                inst_arrays[k] = v
+        return TPInstrument(list(inst_arrays.items()), instrument.evotype, instrument.state_space)
+
+    if to_type in ("full", "static", "static unitary"):
+        from ..operations import convert as _op_convert
+        ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {}
+        members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure))
+                    for k, g in instrument.items()]
+        return Instrument(members, instrument.evotype, instrument.state_space)
+
+    if to_type == 'CPTPLND':
+        op_arrays = {k: v.to_dense('HilbertSchmidt') for (k,v) in instrument.items()}
+        inst = cptp_instrument(op_arrays, basis)
+        return inst
+
+    raise ValueError("Cannot convert an instrument to type %s" % to_type)
+
+
+def cptp_instrument(op_arrays: dict[str, _np.ndarray], basis: Basis,
+                    error_tol: float = 1e-6, trunc_tol: float = 1e-7,
+                    povm_errormap=None) -> Instrument:
+    """
+    Construct a CPTPLND-parameterized instrument from a dictionary of dense superoperator arrays.
+
+    Each operator in `op_arrays` must be a completely positive trace-reducing (CPTR) map.
+    The function Kraus-decomposes each operator, polar-decomposes each Kraus operator into
+    a unitary and a positive semidefinite part, then assembles the result using
+    CPTPLND-parameterized unitaries and a shared `ComposedPOVM` that enforces the
+    completeness (trace-preservation) constraint across all instrument outcomes.
+
+    The construction follows these steps:
+
+    1. For each outcome label, compute a minimal Kraus decomposition and polar-decompose
+       each Kraus operator ``K = u p`` (``p`` PSD, ``u`` unitary).
+    2. Represent each Kraus term as the composition of a ``RootConjOperator`` (encoding
+       ``ρ ↦ p² ρ p²``, parameterized by the POVM effect for ``p²``) and a
+       CPTPLND-parameterized unitary channel.
+    3. Collect all POVM effects into a shared ``ComposedPOVM`` whose error map is
+       constrained to be CPTPLND.  This enforces that the effects sum to ≤ I, which
+       ensures the full instrument is CPTR.
+    4. If an outcome requires more than one Kraus term, wrap the summands in a
+       ``SummedOperator``; otherwise use the single term directly.
+
+    Parameters
+    ----------
+    op_arrays : dict of str → numpy.ndarray
+        Mapping from outcome label to dense superoperator matrix (in the given basis).
+        Each matrix must represent a CPTR map; the function raises or warns if the
+        Kraus decomposition is inconsistent with this.
+    basis : Basis or str
+        The operator basis in which the superoperators are expressed (e.g. ``'pp'``).
+    error_tol : float, optional
+        Tolerance for eigenvalue errors in the minimal Kraus decomposition.  Eigenvalues
+        below ``-error_tol`` cause an error to be raised.  Default ``1e-6``.
+    trunc_tol : float, optional
+        Truncation tolerance for the minimal Kraus decomposition.  Eigenvalues between
+        ``-error_tol`` and ``trunc_tol`` are silently set to zero.  Default ``1e-7``.
+    povm_errormap : LinearOperator, optional
+        A CPTPLND-parameterized error map to use as the shared POVM error map.
+        If ``None`` (default), the identity channel is used as the starting point
+        and converted to CPTPLND parameterization.
+
+    Returns
+    -------
+    Instrument
+        An ``Instrument`` whose member operations are CPTR maps sharing a common set
+        of CPTPLND parameters.  The instrument is not itself a ``TPInstrument`` because
+        the completeness constraint is enforced implicitly through the shared POVM.
+    """
+    unitaries = dict()
+    effects   = dict()
+
+    # Step 1. Build CPTPLND-parameterized unitaries and static POVM effects
+    #   from Kraus decompositions of the provided operators.
+    for oplbl, op in op_arrays.items():
+        krausops = _ot.minimal_kraus_decomposition(op, basis, error_tol, trunc_tol)
+        cur_effects   = []
+        cur_unitaries = []
+        for K in krausops:
+            # Define F(rho) := K rho K^\\dagger. If we compute a polar decomposition
+            # K = u p, then we can express F(rho) = u (p rho p) u^\\dagger, which is
+            # a composition of a unitary channel and the "RootConjOperator" of p@p.
+            u, p = _la.polar(K, side='right')
+            u_linop = _op.StaticUnitaryOp(u, basis)  # type: ignore
+            u_cptp = _op.convert(u_linop, 'CPTPLND', basis)
+            cur_unitaries.append(u_cptp)
+            E_superket = _bt.stdmx_to_vec(p @ p, basis)
+            E = _pv.StaticPOVMEffect(E_superket)
+            cur_effects.append(E)
+        effects[oplbl]   = cur_effects
+        unitaries[oplbl] = cur_unitaries
+
+    # Step 2. Build a CPTPLND-parameterized POVM from the static POVM effects.
+    # 
+    #   We can't use povms.convert(...) because we might have more
+    #   effects than the Hilbert space dimension.
+    #
+    base_povm_effects = {}
+    for oplbl, cur_effects in effects.items():
+        for i, E in enumerate(cur_effects):
+            elbl = Label((oplbl, i))
+            base_povm_effects[elbl] = E
+    base_povm = _BasePOVM(base_povm_effects, preserve_sum=False)
+    udim = int(_np.round(basis.dim ** 0.5))
+    if povm_errormap is None:
+        I_ideal   = _op.StaticUnitaryOp(_np.eye(udim), basis) # type: ignore
+        I_cptplnd = _op.convert(I_ideal, 'CPTPLND', basis) 
+        povm_errormap = I_cptplnd.factorops[1]
+    povm_cptp = _pv.ComposedPOVM(povm_errormap, base_povm)
+
+    # Step 3. Assemble the CPTPLND-parameterized unitaries and POVM
+    #   effects to represent each operator as a completely-positive
+    #   trace-reducing map.
+    # 
+    inst_ops = dict()
+    for lbl in op_arrays:
+        cur_effects   = effects[lbl]
+        cur_unitaries = unitaries[lbl]
+
+        op_summands = []
+        for i, U_cptplnd in enumerate(cur_unitaries):
+            E_cptp = povm_cptp[Label((lbl, i))]
+            rcop = RootConjOperator(E_cptp, basis)
+            cptr = _op.ComposedOp([rcop, U_cptplnd])
+            op_summands.append(cptr)
+        if len(op_summands) == 1:
+            op = op_summands[0]
+        else:
+            op = SummedOperator(op_summands, basis)
 
-    raise ValueError("Could not convert instrument to to type(s): %s" % str(to_types))
+        inst_ops[lbl] = op
+    inst = Instrument(inst_ops)
+    return inst

pygsti/modelmembers/instruments/instrument.py

@@ -19,6 +19,9 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.tools import matrixtools as _mt
 from pygsti.tools import slicetools as _slct
+from pygsti.tools import basistools as _bt
+from pygsti.tools import optools as _ot
+from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.label import Label as _Label
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 
@@ -378,13 +381,14 @@ def acton(self, state):
 
         staterep = state._rep
         outcome_probs_and_states = _collections.OrderedDict()
+
         for lbl, element in self.items():
             output_rep = element._rep.acton(staterep)
-            output_unnormalized_state = output_rep.to_dense()
-            prob = output_unnormalized_state[0] * state.dim**0.25
-            output_normalized_state = output_unnormalized_state / prob  # so [0]th == 1/state_dim**0.25
-            outcome_probs_and_states[lbl] = (prob, _state.StaticState(output_normalized_state, self.evotype,
-                                                                      self.state_space))
+            unnormalized_state_array = output_rep.to_dense()
+            prob = _ot.superket_trace(unnormalized_state_array, element.basis)
+            output_state_array = unnormalized_state_array / prob
+            output_state = _state.StaticState(output_state_array, self.evotype, self.state_space)
+            outcome_probs_and_states[lbl] = (prob, output_state)
 
         return outcome_probs_and_states
 

pygsti/modelmembers/operations/__init__.py

@@ -39,11 +39,13 @@
 from .stochasticop import StochasticNoiseOp
 from .lindbladcoefficients import LindbladCoefficientBlock as _LindbladCoefficientBlock
 from .affineshiftop import AffineShiftOp
+from .cptrop import RootConjOperator, SummedOperator
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.tools import basistools as _bt
 from pygsti.tools import optools as _ot
 from pygsti import SpaceT
 
+
 def create_from_unitary_mx(unitary_mx, op_type, basis='pp', stdname=None, evotype='default', state_space=None):
     """ TODO: docstring - note that op_type can be a list/tuple of types in order of precedence """
     op_type_preferences = verbose_type_from_op_type(op_type)