Merge pull request #255 from pyGSTio/feature-kraus-api

Feature kraus api

pygsti/algorithms/contract.py

@@ -287,7 +287,7 @@ def _contract_to_cp_direct(model, verbosity, tp_also=False, maxiter=100000, tol=
             it += 1
             if it > maxiter: break
 
-        mdl.operations[opLabel] = _op.FullArbitraryOp(new_op, mdl.evotype, mdl.state_space)
+        mdl.operations[opLabel] = _op.FullArbitraryOp(new_op, evotype=mdl.evotype, state_space=mdl.state_space)
 
         if it > maxiter:
             printer.warning("Max iterations exceeded in contract_to_cp_direct")

pygsti/evotypes/basereps.py

@@ -14,6 +14,12 @@
 
 import numpy as _np
 
+
+class POVMRep:
+    """ The base class for all POVM representation classes """
+    pass
+
+
 try:
     from .basereps_cython import OpRep, StateRep, EffectRep, TermRep, PolynomialRep
 except ImportError:

pygsti/evotypes/chp/__init__.py

@@ -10,6 +10,18 @@
 # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
 #***************************************************************************************************
 
+chpexe = None
+
+
+def chpexe_path():
+    from pathlib import Path as _Path
+    if chpexe is None:
+        raise ValueError(("To use 'chp' evotype, please set `pygsti.evotypes.chp.chpexe`"
+                          "to the path to your chp executable."))
+    return _Path(chpexe)
+
+
+from .povmreps import *
 from .effectreps import *
 from .opreps import *
 from .statereps import *

pygsti/evotypes/chp/effectreps.py

@@ -1,5 +1,5 @@
 """
-POVM effect representation classes for the `stabilizer_slow` evolution type.
+POVM effect representation classes for the `chp` evolution type.
 """
 #***************************************************************************************************
 # Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).

pygsti/evotypes/chp/opreps.py

@@ -13,7 +13,7 @@
 import numpy as _np
 from numpy.random import RandomState as _RandomState
 
-from .statereps import _update_chp_op
+from .statereps import _update_chp_op, StateRep as _StateRep
 from .. import basereps as _basereps
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from ...tools import internalgates as _itgs
@@ -34,13 +34,33 @@ def __init__(self, chp_ops, state_space):
     def num_qubits(self):
         return self.state_space.num_qubits
 
-    def chp_ops(self, seed_or_state=None):
+    def acton(self, state):
+        return self.acton_random(state, None)  # default is to ignore rand_state
+
+    def adjoint_acton(self, state):
+        return self.adjoint_acton_random(state, None)  # default is to ignore rand_state
+
+    def acton_random(self, state, rand_state):
+        return _StateRep(state.chp_ops + self._chp_ops(rand_state), state.state_space)
+
+    def adjoint_acton_random(self, state, rand_state):
+        raise NotImplementedError()
+
+    def _chp_ops(self, seed_or_state=None):
         return self.base_chp_ops
 
-    def chp_str(self, seed_or_state=None):
-        op_str = '\n'.join(self.chp_ops(seed_or_state=seed_or_state))
-        if len(op_str) > 0: op_str += '\n'
-        return op_str
+    #def chp_str(self, seed_or_state=None):
+    #    op_str = '\n'.join(self.chp_ops(seed_or_state=seed_or_state))
+    #    if len(op_str) > 0: op_str += '\n'
+    #    return op_str
+
+    def to_dense(self, on_space):
+        try:
+            str_ops = str(self._chp_ops())
+        except Exception:
+            str_ops = self.__class__.__name__  # fall back to class name
+        raise NotImplementedError("CHP op '%s' cannot convert themselves to dense %s-space matrices (yet)" %
+                                  (str_ops, on_space))
 
 
 class OpRepClifford(OpRep):
@@ -83,11 +103,10 @@ def __init__(self, factor_op_reps, state_space):
     def reinit_factor_op_reps(self, factor_op_reps):
         self.factors_reps = factor_op_reps
 
-    def chp_ops(self, seed_or_state=None):
+    def _chp_ops(self, seed_or_state=None):
         ops = []
         for factor in self.factor_reps:
-            ops.extend(factor.chp_ops(seed_or_state=seed_or_state))
-
+            ops.extend(factor._chp_ops(seed_or_state=seed_or_state))
         return ops
 
 
@@ -113,45 +132,79 @@ def __init__(self, state_space, target_labels, embedded_rep):
         self.embedded_to_local_qubit_indices = {str(i): str(j) for i, j in enumerate(qubit_indices)}
 
         # TODO: This doesn't work as nicely for the stochastic op, where chp_ops can be reset between chp_str calls
-        chp_ops = [_update_chp_op(op, self.embedded_to_local_qubit_indices) for op in self.embedded_rep.chp_ops()]
+        chp_ops = [_update_chp_op(op, self.embedded_to_local_qubit_indices) for op in self.embedded_rep._chp_ops()]
         super(OpRepEmbedded, self).__init__(chp_ops, state_space)
 
-    def chp_ops(self, seed_or_state=None):
-        return [_update_chp_op(op, self.embedded_to_local_qubit_indices) for op in self.embedded_rep.chp_ops(seed_or_state=seed_or_state)]
+    def _chp_ops(self, seed_or_state=None):
+        return [_update_chp_op(op, self.embedded_to_local_qubit_indices) for op in self.embedded_rep._chp_ops(seed_or_state=seed_or_state)]
 
 
 class OpRepRepeated(OpRep):
     def __init__(self, rep_to_repeat, num_repetitions, state_space):
         state_space = _StateSpace.cast(state_space)
         self.repeated_rep = rep_to_repeat
         self.num_repetitions = num_repetitions
-        super(OpRepRepeated, self).__init__(self.repeated_rep.chp_ops() * self.num_repetitions, state_space)
-
+        super(OpRepRepeated, self).__init__(self.repeated_rep._chp_ops() * self.num_repetitions, state_space)
 
-class OpRepStochastic(OpRep):
 
-    def __init__(self, basis, rate_poly_dicts, initial_rates, seed_or_state, state_space):
+class OpRepRandomUnitary(OpRep):
 
+    def __init__(self, basis, unitary_rates, unitary_reps, seed_or_state, state_space):
         self.basis = basis
         assert (basis.name in ['pp', 'PP']), "Only Pauli basis is allowed for 'chp' evotype"
 
+        state_space = _StateSpace.cast(state_space)
+        nqubits = state_space.num_qubits
+        assert(self.basis.dim == 4**nqubits), "Must have an integral number of qubits"
+
+        self.unitary_reps = unitary_reps  # superop reps in this evotype (must be reps of *this* evotype)
+        if unitary_rates is None:
+            self.unitary_rates = _np.ones(len(self.unitary_reps), 'd')
+        else:
+            assert(len(unitary_rates) == len(self.unitary_reps))
+            self.unitary_rates = _np.array(unitary_rates, 'd')
+
         if isinstance(seed_or_state, _RandomState):
-            self.rand_state = seed_or_state
+            self._rand_state = seed_or_state
+        elif isinstance(seed_or_state, int):
+            self._rand_state = _RandomState(seed_or_state)
         else:
-            self.rand_state = _RandomState(seed_or_state)
+            self._rand_state = None
 
-        #TODO: need to fix this: `basis` above functions as basis to make superoperators out of, but here we have
-        # a CHP stochastic op which is given a basis for the space - e.g. a dim=2 vector space for 1 qubit, so
-        # we need to distinguish/specify the basis better for this... and what about rate_poly_dicts (see svterm)
-        nqubits = state_space.num_qubits
-        assert(self.basis.dim == 4**nqubits), "Must have an integral number of qubits"
+        super(OpRepRandomUnitary, self).__init__([], state_space)  # don't store any chp_ops in base
 
-        std_chp_ops = _itgs.standard_gatenames_chp_conversions()
+    def update_unitary_rates(self, rates):
+        self.unitary_rates[:] = rates
+
+    def _chp_ops(self, seed_or_state=None):
+        # Optionally override RNG for this call
+        if seed_or_state is not None:
+            if isinstance(seed_or_state, _np.random.RandomState):
+                rand_state = seed_or_state
+            else:
+                rand_state = _np.random.RandomState(seed_or_state)
+        else:
+            rand_state = self._rand_state
+        assert(rand_state is not None), "Forward sim must provide a random state to use CHP and OpRepRandomUnitary"
+
+        rates = self.unitary_rates
+        index = rand_state.choice(len(self.unitary_rates), p=rates)
+        rep = self.unitary_reps[index]
+        return rep._chp_ops()
+
+
+class OpRepStochastic(OpRepRandomUnitary):
+
+    def __init__(self, stochastic_basis, basis, initial_rates, seed_or_state, state_space):
+
+        self.rates = initial_rates
+        self.stochastic_basis = stochastic_basis
+        assert(self.stochastic_basis.first_element_is_identity)
 
-        # For CHP, need to make a Composed + EmbeddedOp for the super operators
-        # For lower overhead, make this directly using the rep instead of with objects
-        self.stochastic_superop_reps = []
-        for label in self.basis.labels[1:]:
+        rates = [1 - sum(initial_rates)] + list(initial_rates)
+        reps = [OpRep([], state_space)]  # start with identity
+        std_chp_ops = _itgs.standard_gatenames_chp_conversions()
+        for label in self.stochastic_basis.labels[1:]:
             combined_chp_ops = []
 
             for i, pauli in enumerate(label):
@@ -160,31 +213,67 @@ def __init__(self, basis, rate_poly_dicts, initial_rates, seed_or_state, state_s
                 chp_op_targeted = [op.replace('0', str(i)) for op in chp_op]
                 combined_chp_ops.extend(chp_op_targeted)
 
-            sub_rep = OpRep(combined_chp_ops, state_space)
-            self.stochastic_superop_reps.append(sub_rep)
-        self.rates = initial_rates
-        super(OpRepStochastic, self).__init__([], state_space)  # don't store any chp_ops in base
+            reps.append(OpRep(combined_chp_ops, state_space))
+
+        super(OpRepStochastic, self).__init__(basis, _np.array(rates, 'd'), reps, seed_or_state, state_space)
+
+        #OLD
+        #self.basis = basis
+        #assert (basis.name == 'pp'), "Only Pauli basis is allowed for 'chp' evotype"
+        #
+        #if isinstance(seed_or_state, _RandomState):
+        #    self.rand_state = seed_or_state
+        #else:
+        #    self.rand_state = _RandomState(seed_or_state)
+        #
+        ##TODO: need to fix this: `basis` above functions as basis to make superoperators out of, but here we have
+        ## a CHP stochastic op which is given a basis for the space - e.g. a dim=2 vector space for 1 qubit, so
+        ## we need to distinguish/specify the basis better for this... and what about rate_poly_dicts (see svterm)
+        #nqubits = state_space.num_qubits
+        #assert(self.basis.dim == 4**nqubits), "Must have an integral number of qubits"
+        #
+        #std_chp_ops = _itgs.standard_gatenames_chp_conversions()
+        #
+        ## For CHP, need to make a Composed + EmbeddedOp for the super operators
+        ## For lower overhead, make this directly using the rep instead of with objects
+        #self.stochastic_superop_reps = []
+        #for label in self.basis.labels[1:]:
+        #    combined_chp_ops = []
+        #
+        #    for i, pauli in enumerate(label):
+        #        name = 'Gi' if pauli == "I" else 'G%spi' % pauli.lower()
+        #        chp_op = std_chp_ops[name]
+        #        chp_op_targeted = [op.replace('0', str(i)) for op in chp_op]
+        #        combined_chp_ops.extend(chp_op_targeted)
+        #
+        #    sub_rep = OpRep(combined_chp_ops, state_space)
+        #    self.stochastic_superop_reps.append(sub_rep)
+        #self.rates = initial_rates
+        #super(OpRepStochastic, self).__init__([], state_space)  # don't store any chp_ops in base
 
     def update_rates(self, rates):
+        unitary_rates = [1 - sum(rates)] + list(rates)
         self.rates[:] = rates
-
-    def chp_ops(self, seed_or_state=None):
-        # Optionally override RNG for this call
-        if seed_or_state is not None:
-            if isinstance(seed_or_state, _np.random.RandomState):
-                rand_state = seed_or_state
-            else:
-                rand_state = _np.random.RandomState(seed_or_state)
-        else:
-            rand_state = self.rand_state
-
-        rates = self.rates
-        all_rates = [*rates, 1.0 - sum(rates)]  # Include identity so that probabilities are 1
-        index = rand_state.choice(self.basis.size, p=all_rates)
-
-        # If final entry, no operation selected
-        if index == self.basis.size - 1:
-            return []
-
-        rep = self.stochastic_superop_reps[index]
-        return rep.chp_ops()
+        self.update_unitary_rates(unitary_rates)
+
+    #TODO REMOVE - covered by OpRepRandomUnitary
+    #def chp_ops(self, seed_or_state=None):
+    #    # Optionally override RNG for this call
+    #    if seed_or_state is not None:
+    #        if isinstance(seed_or_state, _np.random.RandomState):
+    #            rand_state = seed_or_state
+    #        else:
+    #            rand_state = _np.random.RandomState(seed_or_state)
+    #    else:
+    #        rand_state = self.rand_state
+    #
+    #    rates = self.rates
+    #    all_rates = [*rates, 1.0 - sum(rates)]  # Include identity so that probabilities are 1
+    #    index = rand_state.choice(self.basis.size, p=all_rates)
+    #
+    #    # If final entry, no operation selected
+    #    if index == self.basis.size - 1:
+    #        return []
+    #
+    #    rep = self.stochastic_superop_reps[index]
+    #    return rep._chp_ops()