WIP: CPTP instruments

pygsti/modelmembers/instruments/__init__.py

@@ -13,8 +13,17 @@
 from .instrument import Instrument
 from .tpinstrument import TPInstrument
 from .tpinstrumentop import TPInstrumentOp
-
-from pygsti.tools import optools as _ot
+from .cptpinstrument import RootConjOperator, SummedOperator
+
+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,24 +46,19 @@ def instrument_type_from_op_type(op_type):
 
     # Limited set (only matching what is in convert)
     instr_conversion = {
-        'auto': 'full',
         'static unitary': 'static unitary',
         'static clifford': 'static clifford',
         'static': 'static',
         'full': 'full',
         'full TP': 'full TP',
         'full CPTP': 'full CPTP',
         'full unitary': 'full unitary',
+        '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
@@ -110,26 +114,103 @@ 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:
+    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/cptpinstrument.py

@@ -0,0 +1,150 @@
+"""
+Defines the CPTPInstrument class
+"""
+#***************************************************************************************************
+# Copyright 2015, 2019, 2025 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
+# in this software.
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
+# in compliance with the License.  You may obtain a copy of the License at
+# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
+#***************************************************************************************************
+
+import numpy as _np
+from pygsti.pgtypes import SpaceT
+from pygsti.baseobjs.basis import Basis
+from pygsti.modelmembers.povms.effect import POVMEffect
+from pygsti.modelmembers.operations import LinearOperator
+from pygsti.tools import optools as _ot
+
+from typing import Union
+BasisLike = Union[Basis, str]
+
+
+class RootConjOperator(LinearOperator):
+    """
+    A linear operator parameterized by a matrix E where 0 ≤ E ≤ 1, whose action takes ρ to E½ ρ E½ .
+
+    Every CPTR map can be obtained by pre-composing a CPTP map with this kind of linear operator.
+    """
+
+    EIGTOL_WARNING = 1e-10
+    EIGTOL_ERROR   = 1e-8
+
+    def __init__(self, effect: POVMEffect, basis: BasisLike):
+        self._basis        = Basis.cast(basis, effect.dim)
+        self._effect       = effect
+        self._state_space  = effect.state_space
+        self._evotype      = effect.evotype
+
+        dim = self._state_space.dim
+        self._rep = self._evotype.create_dense_superop_rep( _np.zeros((dim, dim)), self._basis, self._state_space )
+        self._update_rep_base()
+        LinearOperator.__init__(self, self._rep, self._evotype)
+        self.init_gpindices()
+
+    def submembers(self):
+        return [self._effect]
+
+    def _update_rep_base(self):
+        # This function is directly analogous to TPInstrumentOp._construct_matrix.
+        self._rep.base.flags.writeable = True
+        assert(self._rep.base.shape == (self.dim, self.dim))
+        effect_superket = self._effect.to_dense()
+        mx = _ot.rootconj_superop(effect_superket, self._basis)
+        self._rep.base[:] = mx
+        self._rep.base.flags.writeable = False
+        self._rep.base_has_changed()
+        return
+
+    def deriv_wrt_params(self, wrt_filter=None):
+        raise NotImplementedError()
+
+    def has_nonzero_hessian(self):
+        # This is not affine in its parameters.
+        return True
+
+    def from_vector(self, v, close=False, dirty_value=True):
+        for sm, local_inds in zip(self.submembers(), self._submember_rpindices):
+            sm.from_vector(v[local_inds], close, dirty_value)
+        self._update_rep_base()
+        return
+
+    @property
+    def num_params(self):
+        return len(self.gpindices_as_array())
+
+    def to_vector(self):
+        v = _np.empty(self.num_params, 'd')
+        for param_op, local_inds in zip(self.submembers(), self._submember_rpindices):
+            v[local_inds] = param_op.to_vector()
+        return v
+
+    def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray:
+        assert on_space in ('HilbertSchmidt', 'minimal')
+        out = self._rep.base.copy()
+        out.flags.writeable = True
+        return out
+
+
+class SummedOperator(LinearOperator):
+
+
+    def __init__(self, operators, basis: BasisLike):
+        op = operators[0]
+        self._basis        = Basis.cast(basis, op.dim)
+        self._operators    = operators
+        self._state_space  = op.state_space
+        self._evotype      = op.evotype
+        self._subreps      = [op._rep for op in self._operators]
+        self._rep = self._evotype.create_sum_rep( self._subreps, self._state_space )
+        LinearOperator.__init__(self, self._rep, self._evotype)
+        self.init_gpindices()
+        # NOTE: This class doesn't have a function analogous to _update_rep_base
+        # that we use in RootConjOperator. We can get away with not having such
+        # a function because it's the responsibility of op.from_vector(...)
+        # to update op's attached OpRep.
+        return
+
+    def submembers(self):
+        out = []
+        hit = set()
+        for op in self._operators:
+            temp = op.submembers()
+            for sm in temp:
+                if id(temp) not in hit:
+                    hit.add(id(temp))
+                    out.append(sm)
+        return out
+
+    def deriv_wrt_params(self, wrt_filter=None):
+        raise NotImplementedError()
+
+    def has_nonzero_hessian(self):
+        return any(op.has_nonzero_hession() for op in self._operators)
+
+    def from_vector(self, v, close=False, dirty_value=True):
+        for sm, local_inds in zip(self.submembers(), self._submember_rpindices):
+            sm.from_vector(v[local_inds], close, dirty_value)
+        return
+
+    @property
+    def num_params(self):
+        return len(self.gpindices_as_array())
+
+    def to_vector(self):
+        v = _np.empty(self.num_params, 'd')
+        for param_op, local_inds in zip(self.submembers(), self._submember_rpindices):
+            v[local_inds] = param_op.to_vector()
+        return v
+
+    def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray:
+        assert on_space in ('HilbertSchmidt', 'minimal')
+        on_space = 'HilbertSchmidt'
+        out = self._operators[0].to_dense(on_space)
+        if not out.flags.writeable:
+            out = out.copy()
+        for op in self._operators[1:]:
+            temp = op.to_dense(on_space)
+            out += temp
+        return out

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