diff --git a/pygsti/__init__.py b/pygsti/__init__.py
index 90a44370a..09bb0c30f 100644
--- a/pygsti/__init__.py
+++ b/pygsti/__init__.py
@@ -7,6 +7,9 @@
 # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
 #***************************************************************************************************
 """ A Python implementation of LinearOperator Set Tomography """
+from .pgtypes import (
+    SpaceT
+)
 
 from . import baseobjs
 from . import algorithms as alg
@@ -19,6 +22,7 @@
 from . import protocols
 from . import report as rpt
 from . import serialization
+from . import enums
 
 # Import the most important/useful routines of each module/sub-package
 # into the package namespace
@@ -30,6 +34,7 @@
 from pygsti.tools.gatetools import *  # *_qubit_gate fns
 from .drivers import *
 from .tools import *
+
 # NUMPY BUG FIX (imported from tools)
 from pygsti.baseobjs._compatibility import _numpy14einsumfix
 
diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py
index 9b3eead56..941e72845 100644
--- a/pygsti/algorithms/gaugeopt.py
+++ b/pygsti/algorithms/gaugeopt.py
@@ -517,8 +517,8 @@ def _jacobian_fn(gauge_group_el):
                 # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS)
                 #   Note: (S_inv * G * S) is G' (transformed G)
                 wt = item_weights.get(lbl, opWeight)
-                left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense(on_space='minimal'))  # shape (n,d1,d2)
-                right = _np.swapaxes(_np.dot(G.to_dense(on_space='minimal'), dS), 0, 1)  # shape (d1,n,d2) -> (n,d1,d2)
+                left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal'))  # shape (n,d1,d2)
+                right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1)  # shape (d1,n,d2) -> (n,d1,d2)
                 result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2)  # shape (d1, d2, n)
                 result = result.reshape((d**2, n))  # must copy b/c non-contiguous
                 my_jacMx[start:start + d**2] = wt * result
@@ -530,8 +530,8 @@ def _jacobian_fn(gauge_group_el):
                 wt = item_weights.get(ilbl, opWeight)
                 for lbl, G in Inst.items():
                     # same calculation as for operation terms
-                    left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense(on_space='minimal'))  # (n,d1,d2)
-                    right = _np.swapaxes(_np.dot(G.to_dense(on_space='minimal'), dS), 0, 1)  # (d1,n,d2) -> (n,d1,d2)
+                    left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal'))  # (n,d1,d2)
+                    right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1)  # (d1,n,d2) -> (n,d1,d2)
                     result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2)  # shape (d1, d2, n)
                     result = result.reshape((d**2, n))  # must copy b/c non-contiguous
                     my_jacMx[start:start + d**2] = wt * result
@@ -544,7 +544,7 @@ def _jacobian_fn(gauge_group_el):
                 #   Note: (S_inv * rho) is transformed rho
                 wt = item_weights.get(lbl, spamWeight)
                 Sinv_dS = _np.dot(S_inv, dS)  # shape (d1,n,d2)
-                result = -1 * _np.dot(Sinv_dS, rho.to_dense(on_space='minimal'))  # shape (d,n)
+                result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal'))  # shape (d,n)
                 my_jacMx[start:start + d] = wt * result
                 start += d
 
@@ -554,7 +554,7 @@ def _jacobian_fn(gauge_group_el):
                 for lbl, E in povm.items():
                     # d(ET_term) = E.T * dS
                     wt = item_weights.get(povmlbl + "_" + lbl, spamWeight)
-                    result = _np.dot(E.to_dense(on_space='minimal')[None, :], dS).T  # shape (1,n,d2).T => (d2,n,1)
+                    result = _np.dot(E.to_dense('minimal')[None, :], dS).T  # shape (1,n,d2).T => (d2,n,1)
                     my_jacMx[start:start + d] = wt * result.squeeze(2)  # (d2,n)
                     start += d
 
@@ -851,7 +851,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl_pre, mdl_post,
 
         #get sgn(denMx) == d(|denMx|_Tr)/d(denMx) in std basis
         # dmDim = denMx.shape[0]
-        denMx = _tools.vec_to_stdmx(prepvec.to_dense(on_space='minimal')[:, None], op_basis)
+        denMx = _tools.vec_to_stdmx(prepvec.to_dense('minimal')[:, None], op_basis)
         assert(_np.linalg.norm(denMx - denMx.T.conjugate()) < 1e-4), \
             "denMx should be Hermitian!"
 
@@ -865,7 +865,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl_pre, mdl_post,
         # get d(prepvec')/dp = d(S_inv * prepvec)/dp in op_basis [shape == (n,dim)]
         #                    = (-S_inv*dS*S_inv) * prepvec = -S_inv*dS * prepvec'
         Sinv_dS = _np.dot(S_inv, dS)  # shape (d1,n,d2)
-        dVdp = -1 * _np.dot(Sinv_dS, prepvec.to_dense(on_space='minimal')[:, None]).squeeze(2)  # shape (d,n,1) => (d,n)
+        dVdp = -1 * _np.dot(Sinv_dS, prepvec.to_dense('minimal')[:, None]).squeeze(2)  # shape (d,n,1) => (d,n)
         assert(dVdp.shape == (d, n))
 
         # denMx = sum( spamvec[i] * Bmx[i] )
@@ -890,7 +890,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl_pre, mdl_post,
         for lbl, effectvec in povm.items():
 
             #get sgn(EMx) == d(|EMx|_Tr)/d(EMx) in std basis
-            EMx = _tools.vec_to_stdmx(effectvec.to_dense(on_space='minimal')[:, None], op_basis)
+            EMx = _tools.vec_to_stdmx(effectvec.to_dense('minimal')[:, None], op_basis)
             # dmDim = EMx.shape[0]
             assert(_np.linalg.norm(EMx - EMx.T.conjugate()) < 1e-4), \
                 "denMx should be Hermitian!"
@@ -905,7 +905,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl_pre, mdl_post,
             # get d(effectvec')/dp = [d(effectvec.T * S)/dp].T in op_basis [shape == (n,dim)]
             #                      = [effectvec.T * dS].T
             #  OR = dS.T * effectvec
-            pre_effectvec = mdl_pre.povms[povmlbl][lbl].to_dense(on_space='minimal')[:, None]
+            pre_effectvec = mdl_pre.povms[povmlbl][lbl].to_dense('minimal')[:, None]
             dVdp = _np.dot(pre_effectvec.T, dS).squeeze(0).T
             # shape = (1,d) * (n, d1,d2) = (1,n,d2) => (n,d2) => (d2,n)
             assert(dVdp.shape == (d, n))
diff --git a/pygsti/circuits/cloudcircuitconstruction.py b/pygsti/circuits/cloudcircuitconstruction.py
index c35b3c608..a7ee63a82 100644
--- a/pygsti/circuits/cloudcircuitconstruction.py
+++ b/pygsti/circuits/cloudcircuitconstruction.py
@@ -27,11 +27,8 @@
 from pygsti.circuits.circuitstructure import GermFiducialPairPlaquette as _GermFiducialPairPlaquette, \
     PlaquetteGridCircuitStructure as _PlaquetteGridCircuitStructure
 
-from pygsti.tools import basistools as _bt
-from pygsti.tools import internalgates as _itgs
 from pygsti.tools import listtools as _lt
 from pygsti.tools import mpitools as _mpit
-from pygsti.tools import optools as _ot
 from pygsti.tools import slicetools as _slct
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
 
diff --git a/pygsti/enums/__init__.py b/pygsti/enums/__init__.py
new file mode 100644
index 000000000..d7558f20b
--- /dev/null
+++ b/pygsti/enums/__init__.py
@@ -0,0 +1,6 @@
+"""
+These enums are intended to make it simplier to call functions with fixed strings expected as input.
+"""
+
+
+from .convertspaceenum import SpaceConversionType
\ No newline at end of file
diff --git a/pygsti/enums/convertspaceenum.py b/pygsti/enums/convertspaceenum.py
new file mode 100644
index 000000000..f774351f9
--- /dev/null
+++ b/pygsti/enums/convertspaceenum.py
@@ -0,0 +1,7 @@
+from enum import Enum
+
+
+class SpaceConversionType(Enum):
+    Minimal = 'minimal'
+    Hilbert = 'Hilbert'
+    HilbertSchmidt = 'HilbertSchmidt'
\ No newline at end of file
diff --git a/pygsti/evotypes/chp/opreps.py b/pygsti/evotypes/chp/opreps.py
index f042b2244..c21f728a8 100644
--- a/pygsti/evotypes/chp/opreps.py
+++ b/pygsti/evotypes/chp/opreps.py
@@ -17,6 +17,7 @@
 from .. import basereps as _basereps
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from ...tools import internalgates as _itgs
+from pygsti import SpaceT
 
 
 class OpRep(_basereps.OpRep):
@@ -49,7 +50,7 @@ def adjoint_acton_random(self, state, rand_state):
     def _chp_ops(self, seed_or_state=None):
         return self.base_chp_ops
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         try:
             str_ops = str(self._chp_ops())
         except Exception:
diff --git a/pygsti/evotypes/densitymx_slow/effectreps.py b/pygsti/evotypes/densitymx_slow/effectreps.py
index 62c2e358c..1ca8490a2 100644
--- a/pygsti/evotypes/densitymx_slow/effectreps.py
+++ b/pygsti/evotypes/densitymx_slow/effectreps.py
@@ -15,7 +15,7 @@
 # import functools as _functools
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from ...tools import matrixtools as _mt
-
+from pygsti import SpaceT
 
 class EffectRep:
     """Any representation of an "effect" in the sense of a POVM."""
@@ -44,7 +44,7 @@ def probability(self, state):
         # can assume state is a StateRep and self.state_rep is
         return _np.dot(self.state_rep.data, state.data)  # not vdot b/c *real* data
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         return self.state_rep.to_dense(on_space)
 
 
@@ -79,7 +79,7 @@ def probability(self, state):
         Edense = self.to_dense('HilbertSchmidt', scratch)
         return _np.dot(Edense, state.data)  # not vdot b/c data is *real*
 
-    def to_dense(self, on_space, outvec=None):
+    def to_dense(self, on_space: SpaceT, outvec=None):
         if on_space not in ('minimal', 'HilbertSchmidt'):
             raise ValueError("'densitymx' evotype cannot produce Hilbert-space ops!")
         return _mt.zvals_int64_to_dense(self.zvals_int, self.nfactors, outvec, False, self.abs_elval)
@@ -105,7 +105,7 @@ def __init__(self, povm_factors, effect_labels, state_space):
         super(EffectRepTensorProduct, self).__init__(state_space)
         self.factor_effects_have_changed()
 
-    def to_dense(self, on_space, outvec=None):
+    def to_dense(self, on_space: SpaceT, outvec=None):
 
         if on_space not in ('minimal', 'HilbertSchmidt'):
             raise ValueError("'densitymx' evotype cannot produce Hilbert-space ops!")
diff --git a/pygsti/evotypes/densitymx_slow/opreps.py b/pygsti/evotypes/densitymx_slow/opreps.py
index af3ffc5b2..0a1d46afd 100644
--- a/pygsti/evotypes/densitymx_slow/opreps.py
+++ b/pygsti/evotypes/densitymx_slow/opreps.py
@@ -23,7 +23,7 @@
 from ...tools import lindbladtools as _lbt
 from ...tools import matrixtools as _mt
 from ...tools import optools as _ot
-
+from pygsti import SpaceT
 
 class OpRep:
     """
@@ -51,12 +51,12 @@ def aslinearoperator(self):
         def mv(v):
             if v.ndim == 2 and v.shape[1] == 1: v = v[:, 0]
             in_state = _StateRepDense(_np.ascontiguousarray(v, 'd'), self.state_space, None)
-            return self.acton(in_state).to_dense(on_space='HilbertSchmidt')
+            return self.acton(in_state).to_dense(on_space="HilbertSchmidt")
 
         def rmv(v):
             if v.ndim == 2 and v.shape[1] == 1: v = v[:, 0]
             in_state = _StateRepDense(_np.ascontiguousarray(v, 'd'), self.state_space, None)
-            return self.adjoint_acton(in_state).to_dense(on_space='HilbertSchmidt')
+            return self.adjoint_acton(in_state).to_dense(on_space="HilbertSchmidt")
         return LinearOperator((self.dim, self.dim), matvec=mv, rmatvec=rmv)  # transpose, adjoint, dot, matmat?
 
 
@@ -80,7 +80,7 @@ def __init__(self, mx, basis, state_space):
     def base_has_changed(self):
         pass
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         if on_space not in ('minimal', 'HilbertSchmidt'):
             raise ValueError("'densitymx_slow' evotype cannot produce Hilbert-space ops!")
         return self.base
@@ -115,7 +115,7 @@ def __init__(self, mx, basis, state_space):
     def base_has_changed(self):
         self.superop_base[:, :] = _ot.unitary_to_superop(self.base, self.basis)
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         if on_space in ('minimal', 'HilbertSchmidt'):
             return self.to_dense_superop()
         else:  # 'Hilbert'
@@ -167,7 +167,7 @@ def adjoint_acton(self, state):
         Aadj = self.A.conjugate(copy=True).transpose()
         return _StateRepDense(Aadj.dot(state.data), state.state_space, None)
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         if on_space not in ('minimal', 'HilbertSchmidt'):
             raise ValueError("'densitymx_slow' evotype cannot produce Hilbert-space ops!")
         return self.A.toarray()
@@ -212,7 +212,7 @@ def __str__(self):
     def copy(self):
         return OpRepKraus(self.basis, list(self.kraus_reps), self.state_space)
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         assert(on_space in ('minimal', 'HilbertSchmidt')), \
             'Can only compute OpRepKraus.to_dense on HilbertSchmidt space!'
         return sum([rep.to_dense(on_space) for rep in self.kraus_reps])
@@ -245,7 +245,7 @@ def copy(self):
     def update_unitary_rates(self, rates):
         self.unitary_rates[:] = rates
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         assert(on_space in ('minimal', 'HilbertSchmidt'))  # below code only works in this case
         return sum([rate * rep.to_dense(on_space) for rate, rep in zip(self.unitary_rates, self.unitary_reps)])
 
diff --git a/pygsti/evotypes/densitymx_slow/statereps.py b/pygsti/evotypes/densitymx_slow/statereps.py
index 69d321d1a..cb5e73310 100644
--- a/pygsti/evotypes/densitymx_slow/statereps.py
+++ b/pygsti/evotypes/densitymx_slow/statereps.py
@@ -17,7 +17,7 @@
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from ...tools import basistools as _bt
 from ...tools import optools as _ot
-
+from pygsti import SpaceT
 try:
     from ...tools import fastcalc as _fastcalc
 except ImportError:
@@ -51,7 +51,7 @@ def actionable_staterep(self):
         # a probability/amplitude by POVM effect reps.
         return self  # for most classes, the rep itself is actionable
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         if on_space not in ('minimal', 'HilbertSchmidt'):
             raise ValueError("'densitymx' evotype cannot produce Hilbert-space ops!")
         return self.data
@@ -136,7 +136,7 @@ class StateRepComposed(StateRep):
     def __init__(self, state_rep, op_rep, state_space):
         self.state_rep = state_rep
         self.op_rep = op_rep
-        super(StateRepComposed, self).__init__(state_rep.to_dense('HilbertSchmidt'), state_space)
+        super(StateRepComposed, self).__init__(state_rep.to_dense("HilbertSchmidt"), state_space)
         self.reps_have_changed()
 
     def reps_have_changed(self):
@@ -158,9 +158,9 @@ def reps_have_changed(self):
         if len(self.factor_reps) == 0:
             vec = _np.empty(0, 'd')
         else:
-            vec = self.factor_reps[0].to_dense('HilbertSchmidt')
+            vec = self.factor_reps[0].to_dense("HilbertSchmidt")
             for i in range(1, len(self.factor_reps)):
-                vec = _np.kron(vec, self.factor_reps[i].to_dense('HilbertSchmidt'))
+                vec = _np.kron(vec, self.factor_reps[i].to_dense("HilbertSchmidt"))
         self.data[:] = vec
 
     def __reduce__(self):
diff --git a/pygsti/evotypes/stabilizer_slow/effectreps.py b/pygsti/evotypes/stabilizer_slow/effectreps.py
index 05797173e..997860686 100644
--- a/pygsti/evotypes/stabilizer_slow/effectreps.py
+++ b/pygsti/evotypes/stabilizer_slow/effectreps.py
@@ -13,7 +13,7 @@
 from .. import basereps as _basereps
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from ...tools import matrixtools as _mt
-
+from pygsti import SpaceT
 
 class EffectRep(_basereps.EffectRep):
     def __init__(self, state_space):
@@ -29,7 +29,7 @@ def probability(self, state):
     def amplitude(self, state):
         return state.sframe.extract_amplitude(self.zvals)
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         return _mt.zvals_to_dense(self.zvals, superket=bool(on_space not in ('minimal', 'Hilbert')))
 
 
@@ -46,7 +46,7 @@ def __str__(self):
         s = "Stabilizer effect vector for %d qubits with outcome %s" % (nQubits, str(self.zvals))
         return s
 
-    def to_dense(self, on_space, outvec=None):
+    def to_dense(self, on_space: SpaceT, outvec=None):
         return _mt.zvals_to_dense(self.zvals, superket=bool(on_space not in ('minimal', 'Hilbert')))
 
 
diff --git a/pygsti/evotypes/statevec_slow/effectreps.py b/pygsti/evotypes/statevec_slow/effectreps.py
index 7e5d01436..c5f584057 100644
--- a/pygsti/evotypes/statevec_slow/effectreps.py
+++ b/pygsti/evotypes/statevec_slow/effectreps.py
@@ -13,7 +13,7 @@
 import numpy as _np
 
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
-
+from pygsti import SpaceT
 
 class EffectRep(object):
     def __init__(self, state_space):
@@ -38,7 +38,7 @@ def __init__(self, state_rep):
     def __str__(self):
         return str(self.state_rep.data)
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         return self.state_rep.to_dense(on_space)
 
     def amplitude(self, state):
@@ -71,7 +71,7 @@ def __init__(self, zvals, basis, state_space):
             base //= 2  # or right shift?
         super(EffectRepComputational, self).__init__(state_space)
 
-    def to_dense(self, on_space, outvec, trust_outvec_sparsity=False):
+    def to_dense(self, on_space: SpaceT, outvec, trust_outvec_sparsity=False):
         # when trust_outvec_sparsity is True, assume we only need to fill in the
         # non-zero elements of outvec (i.e. that outvec is already zero wherever
         # this vector is zero).
@@ -84,7 +84,7 @@ def to_dense(self, on_space, outvec, trust_outvec_sparsity=False):
 
     def amplitude(self, state):  # allow scratch to be passed in?
         scratch = _np.empty(self.dim, complex)
-        Edense = self.to_dense('Hilbert', scratch)
+        Edense = self.to_dense("Hilbert", scratch)
         return _np.vdot(Edense, state.data)
 
 
@@ -113,12 +113,12 @@ def __init__(self, povm_factors, effect_labels, state_space):
     def _fill_fast_kron(self):
         """ Fills in self._fast_kron_array based on current self.factors """
         for i, (factor_dim, Elbl) in enumerate(zip(self._fast_kron_factordims, self.effectLbls)):
-            self.kron_array[i][0:factor_dim] = self.povm_factors[i][Elbl].to_dense('Hilbert')
+            self.kron_array[i][0:factor_dim] = self.povm_factors[i][Elbl].to_dense("Hilbert")
 
     def factor_effects_have_changed(self):
         self._fill_fast_kron()  # updates effect reps
 
-    def to_dense(self, on_space, scratch=None):
+    def to_dense(self, on_space: SpaceT, scratch=None):
         #OLD & SLOW:
         #if len(self.factors) == 0: return _np.empty(0, complex)
         #factorPOVMs = self.factors
@@ -165,7 +165,7 @@ def to_dense(self, on_space, scratch=None):
 
     def amplitude(self, state):  # allow scratch to be passed in?
         scratch = _np.empty(self.dim, complex)
-        Edense = self.to_dense('Hilbert', scratch)
+        Edense = self.to_dense("Hilbert", scratch)
         return _np.vdot(Edense, state.data)
 
 
diff --git a/pygsti/evotypes/statevec_slow/opreps.py b/pygsti/evotypes/statevec_slow/opreps.py
index 0f6108ea2..c9fa8b8c0 100644
--- a/pygsti/evotypes/statevec_slow/opreps.py
+++ b/pygsti/evotypes/statevec_slow/opreps.py
@@ -20,10 +20,9 @@
 from .statereps import StateRepDensePure as _StateRepDensePure
 from .. import basereps as _basereps
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
-from ...tools import basistools as _bt
 from ...tools import internalgates as _itgs
 from ...tools import optools as _ot
-
+from pygsti import SpaceT
 
 class OpRep(_basereps.OpRep):
     def __init__(self, state_space):
@@ -49,12 +48,12 @@ def aslinearoperator(self):
         def mv(v):
             if v.ndim == 2 and v.shape[1] == 1: v = v[:, 0]
             in_state = _StateRepDensePure(_np.ascontiguousarray(v, complex), self.state_space, basis=None)
-            return self.acton(in_state).to_dense('Hilbert')
+            return self.acton(in_state).to_dense("Hilbert")
 
         def rmv(v):
             if v.ndim == 2 and v.shape[1] == 1: v = v[:, 0]
             in_state = _StateRepDensePure(_np.ascontiguousarray(v, complex), self.state_space, basis=None)
-            return self.adjoint_acton(in_state).to_dense('Hilbert')
+            return self.adjoint_acton(in_state).to_dense("Hilbert")
         return LinearOperator((self.dim, self.dim), matvec=mv, rmatvec=rmv)  # transpose, adjoint, dot, matmat?
 
     def copy(self):
@@ -74,7 +73,7 @@ def __init__(self, mx, basis, state_space):
     def base_has_changed(self):
         pass
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         if on_space in ('minimal', 'Hilbert'):
             return self.base
         elif on_space == 'HilbertSchmidt':
@@ -429,7 +428,7 @@ def copy(self):
     def update_unitary_rates(self, rates):
         self.unitary_rates[:] = rates
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         assert(on_space == 'HilbertSchmidt')  # below code only works in this case
         return sum([rate * rep.to_dense(on_space) for rate, rep in zip(self.unitary_rates, self.unitary_reps)])
 
diff --git a/pygsti/evotypes/statevec_slow/statereps.py b/pygsti/evotypes/statevec_slow/statereps.py
index c45783b9e..3a6cfd139 100644
--- a/pygsti/evotypes/statevec_slow/statereps.py
+++ b/pygsti/evotypes/statevec_slow/statereps.py
@@ -18,6 +18,7 @@
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from ...tools import basistools as _bt
 from ...tools import optools as _ot
+from pygsti import SpaceT
 
 try:
     from ...tools import fastcalc as _fastcalc
@@ -48,7 +49,7 @@ def actionable_staterep(self):
         # a probability/amplitude by POVM effect reps.
         return self  # for most classes, the rep itself is actionable
 
-    def to_dense(self, on_space):
+    def to_dense(self, on_space: SpaceT):
         if on_space in ('minimal', 'Hilbert'):
             return self.data
         elif on_space == 'HilbertSchmidt':
@@ -107,7 +108,7 @@ def __init__(self, state_rep, op_rep, state_space):
         self.op_rep = op_rep
         if state_space is None:
             state_space = op_rep.state_space if (op_rep is not None) else state_rep.state_space
-        super(StateRepComposed, self).__init__(state_rep.to_dense('Hilbert'), state_space, self.state_rep.basis)
+        super(StateRepComposed, self).__init__(state_rep.to_dense("Hilbert"), state_space, self.state_rep.basis)
         self.reps_have_changed()
 
     def reps_have_changed(self):
@@ -132,7 +133,7 @@ def reps_have_changed(self):
         if len(self.factor_reps) == 0:
             vec = _np.empty(0, complex)
         else:
-            vec = self.factor_reps[0].to_dense('Hilbert')
+            vec = self.factor_reps[0].to_dense("Hilbert")
             for i in range(1, len(self.factors_reps)):
-                vec = _np.kron(vec, self.factor_reps[i].to_dense('Hilbert'))
+                vec = _np.kron(vec, self.factor_reps[i].to_dense("Hilbert"))
         self.base[:] = vec
diff --git a/pygsti/forwardsims/mapforwardsim.py b/pygsti/forwardsims/mapforwardsim.py
index 81ccf917c..b2f7d3cec 100644
--- a/pygsti/forwardsims/mapforwardsim.py
+++ b/pygsti/forwardsims/mapforwardsim.py
@@ -27,6 +27,7 @@
 from pygsti.tools import slicetools as _slct
 from pygsti.tools.matrixtools import _fas
 from pygsti.tools import listtools as _lt
+from pygsti import SpaceT
 from pygsti.circuits import CircuitList as _CircuitList
 
 _dummy_profiler = _DummyProfiler()
@@ -743,7 +744,7 @@ def product(self, circuit, scale=False):
             G = _np.identity(self.model.evotype.minimal_dim(self.model.state_space))
             for lOp in circuit:
                 if lOp not in scaledGatesAndExps:
-                    opmx = self.model.circuit_layer_operator(lOp, 'op').to_dense(on_space='minimal')
+                    opmx = self.model.circuit_layer_operator(lOp, 'op').to_dense("minimal")
                     ng = max(_nla.norm(opmx), 1.0)
                     scaledGatesAndExps[lOp] = (opmx / ng, _np.log(ng))
 
@@ -764,6 +765,6 @@ def product(self, circuit, scale=False):
         else:
             G = _np.identity(self.model.evotype.minimal_dim(self.model.state_space))
             for lOp in circuit:
-                G = _np.dot(self.model.circuit_layer_operator(lOp, 'op').to_dense(on_space='minimal'), G)
+                G = _np.dot(self.model.circuit_layer_operator(lOp, 'op').to_dense("minimal"), G)
                 # above line: LEXICOGRAPHICAL VS MATRIX ORDER
             return G
\ No newline at end of file
diff --git a/pygsti/forwardsims/matrixforwardsim.py b/pygsti/forwardsims/matrixforwardsim.py
index e811daf73..75b745305 100644
--- a/pygsti/forwardsims/matrixforwardsim.py
+++ b/pygsti/forwardsims/matrixforwardsim.py
@@ -29,6 +29,7 @@
 from pygsti.tools import sharedmemtools as _smt
 from pygsti.tools import slicetools as _slct
 from pygsti.tools.matrixtools import _fas
+from pygsti import SpaceT
 from pygsti.tools import listtools as _lt
 from pygsti.circuits import CircuitList as _CircuitList
 
@@ -86,7 +87,7 @@ def product(self, circuit, scale=False):
             G = _np.identity(self.model.evotype.minimal_dim(self.model.state_space))
             for lOp in circuit:
                 if lOp not in scaledGatesAndExps:
-                    opmx = self.model.circuit_layer_operator(lOp, 'op').to_dense(on_space='minimal')
+                    opmx = self.model.circuit_layer_operator(lOp, 'op').to_dense("minimal")
                     ng = max(_nla.norm(opmx), 1.0)
                     scaledGatesAndExps[lOp] = (opmx / ng, _np.log(ng))
 
@@ -107,15 +108,15 @@ def product(self, circuit, scale=False):
         else:
             G = _np.identity(self.model.evotype.minimal_dim(self.model.state_space))
             for lOp in circuit:
-                G = _np.dot(self.model.circuit_layer_operator(lOp, 'op').to_dense(on_space='minimal'), G)
+                G = _np.dot(self.model.circuit_layer_operator(lOp, 'op').to_dense("minimal"), G)
                 # above line: LEXICOGRAPHICAL VS MATRIX ORDER
             return G
 
     def _rho_es_from_spam_tuples(self, rholabel, elabels):
         # This calculator uses the convention that rho has shape (N,1)
-        rho = self.model.circuit_layer_operator(rholabel, 'prep').to_dense(on_space='minimal')[:, None]
+        rho = self.model.circuit_layer_operator(rholabel, 'prep').to_dense("minimal")[:, None]
         Es = [_np.conjugate(_np.transpose(self.model.circuit_layer_operator(
-              elabel, 'povm').to_dense(on_space='minimal')[:, None]))
+              elabel, 'povm').to_dense("minimal")[:, None]))
               for elabel in elabels]  # [:, None] becuse of convention: E has shape (1,N)
         return rho, Es
 
@@ -341,13 +342,13 @@ def dproduct(self, circuit, flat=False, wrt_filter=None):
         leftProds = []
         G = _np.identity(dim); leftProds.append(G)
         for opLabel in revOpLabelList:
-            G = _np.dot(G, self.model.circuit_layer_operator(opLabel, 'op').to_dense(on_space='minimal'))
+            G = _np.dot(G, self.model.circuit_layer_operator(opLabel, 'op').to_dense("minimal"))
             leftProds.append(G)
 
         rightProdsT = []
         G = _np.identity(dim); rightProdsT.append(_np.transpose(G))
         for opLabel in reversed(revOpLabelList):
-            G = _np.dot(self.model.circuit_layer_operator(opLabel, 'op').to_dense(on_space='minimal'), G)
+            G = _np.dot(self.model.circuit_layer_operator(opLabel, 'op').to_dense("minimal"), G)
             rightProdsT.append(_np.transpose(G))
 
         # Allocate memory for the final result
@@ -468,7 +469,7 @@ def hproduct(self, circuit, flat=False, wrt_filter1=None, wrt_filter2=None):
             prods[(i, i - 1)] = ident  # product of no gates
             G = ident
             for (j, opLabel2) in enumerate(revOpLabelList[i:], start=i):  # loop over "ending" gate (>= starting gate)
-                G = _np.dot(G, self.model.circuit_layer_operator(opLabel2, 'op').to_dense(on_space='minimal'))
+                G = _np.dot(G, self.model.circuit_layer_operator(opLabel2, 'op').to_dense("minimal"))
                 prods[(i, j)] = G
         prods[(len(revOpLabelList), len(revOpLabelList) - 1)] = ident  # product of no gates
 
@@ -749,7 +750,7 @@ def _compute_product_cache(self, layout_atom_tree, resource_alloc):
                     prodCache[iDest] = _np.identity(dim)
                     # Note: scaleCache[i] = 0.0 from initialization
                 else:
-                    gate = self.model.circuit_layer_operator(opLabel, 'op').to_dense(on_space='minimal')
+                    gate = self.model.circuit_layer_operator(opLabel, 'op').to_dense("minimal")
                     nG = max(_nla.norm(gate), 1.0)
                     prodCache[iDest] = gate / nG
                     scaleCache[iDest] = _np.log(nG)
@@ -1205,9 +1206,9 @@ def _scale_exp(self, scale_exps):
     def _rho_e_from_spam_tuple(self, spam_tuple):
         # This calculator uses the convention that rho has shape (N,1)
         rholabel, elabel = spam_tuple
-        rho = self.model.circuit_layer_operator(rholabel, 'prep').to_dense(on_space='minimal')[:, None]
+        rho = self.model.circuit_layer_operator(rholabel, 'prep').to_dense("minimal")[:, None]
         E = _np.conjugate(_np.transpose(self.model.circuit_layer_operator(
-            elabel, 'povm').to_dense(on_space='minimal')[:, None]))
+            elabel, 'povm').to_dense("minimal")[:, None]))
         return rho, E
 
     def _probs_from_rho_e(self, rho, e, gs, scale_vals):
diff --git a/pygsti/io/writers.py b/pygsti/io/writers.py
index d93aab5cb..4ffd0f63c 100644
--- a/pygsti/io/writers.py
+++ b/pygsti/io/writers.py
@@ -23,6 +23,7 @@
 
 # from . import stdinput as _stdinput
 from pygsti import tools as _tools
+from pygsti import SpaceT
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
 from pygsti.modelmembers import instruments as _instrument
 from pygsti.modelmembers import operations as _op
@@ -349,8 +350,8 @@ def writeprop(f, lbl, val):
             elif isinstance(rhoVec, _state.StaticState): typ = "STATIC-PREP"
             #elif isinstance(rhoVec, _state.LindbladSPAMVec):  # TODO - change to ComposedState?
             #    typ = "CPTP-PREP"
-            #    props = [("PureVec", rhoVec.state_vec.to_dense(on_space='HilbertSchmidt')),
-            #             ("ErrgenMx", rhoVec.error_map.to_dense(on_space='HilbertSchmidt'))]
+            #    props = [("PureVec", rhoVec.state_vec.to_dense("HilbertSchmidt")),
+            #             ("ErrgenMx", rhoVec.error_map.to_dense("HilbertSchmidt"))]
             else:
                 _warnings.warn(
                     ("Non-standard prep of type {typ} cannot be described by"
@@ -358,7 +359,7 @@ def writeprop(f, lbl, val):
                      "fully parameterized spam vector").format(typ=str(type(rhoVec))))
                 typ = "PREP"
 
-            if props is None: props = [("LiouvilleVec", rhoVec.to_dense(on_space='HilbertSchmidt'))]
+            if props is None: props = [("LiouvilleVec", rhoVec.to_dense("HilbertSchmidt"))]
             output.write("%s: %s\n" % (typ, prepLabel))
             for lbl, val in props:
                 writeprop(output, lbl, val)
@@ -369,7 +370,7 @@ def writeprop(f, lbl, val):
             elif isinstance(povm, _povm.TPPOVM): povmType = "TP-POVM"
             #elif isinstance(povm, _povm.LindbladPOVM):  # TODO - change to ComposedPOVM?
             #    povmType = "CPTP-POVM"
-            #    props = [("ErrgenMx", povm.error_map.to_dense(on_space='HilbertSchmidt'))]
+            #    props = [("ErrgenMx", povm.error_map.to_dense("HilbertSchmidt"))]
             #    povm_to_write = povm.base_povm
             else:
                 _warnings.warn(
@@ -394,7 +395,7 @@ def writeprop(f, lbl, val):
                          "fully parameterized spam vector").format(typ=str(type(EVec))))
                     typ = "EFFECT"
                 output.write("%s: %s\n" % (typ, ELabel))
-                writeprop(output, "LiouvilleVec", EVec.to_dense(on_space='HilbertSchmidt'))
+                writeprop(output, "LiouvilleVec", EVec.to_dense("HilbertSchmidt"))
 
             output.write("END POVM\n\n")
 
@@ -405,7 +406,7 @@ def writeprop(f, lbl, val):
             elif isinstance(gate, _op.StaticArbitraryOp): typ = "STATIC-GATE"
             elif isinstance(gate, _op.ComposedOp):
                 typ = "COMPOSED-GATE"
-                props = [("%dLiouvilleMx" % i, factor.to_dense(on_space='HilbertSchmidt'))
+                props = [("%dLiouvilleMx" % i, factor.to_dense("HilbertSchmidt"))
                          for i, factor in enumerate(gate.factorops)]
             else:
                 _warnings.warn(
@@ -414,7 +415,7 @@ def writeprop(f, lbl, val):
                      "fully parameterized gate").format(typ=str(type(gate))))
                 typ = "GATE"
 
-            if props is None: props = [("LiouvilleMx", gate.to_dense(on_space='HilbertSchmidt'))]
+            if props is None: props = [("LiouvilleMx", gate.to_dense("HilbertSchmidt"))]
             output.write(typ + ": " + str(label) + '\n')
             for lbl, val in props:
                 writeprop(output, lbl, val)
@@ -441,7 +442,7 @@ def writeprop(f, lbl, val):
                          "fully parameterized gate").format(typ=str(type(gate))))
                     typ = "IGATE"
                 output.write(typ + ": " + str(label) + '\n')
-                writeprop(output, "LiouvilleMx", gate.to_dense(on_space='HilbertSchmidt'))
+                writeprop(output, "LiouvilleMx", gate.to_dense("HilbertSchmidt"))
             output.write("END Instrument\n\n")
 
         if model.state_space is not None:
diff --git a/pygsti/layouts/copalayout.py b/pygsti/layouts/copalayout.py
index 38b647f5a..f2ea0d8e6 100644
--- a/pygsti/layouts/copalayout.py
+++ b/pygsti/layouts/copalayout.py
@@ -87,10 +87,6 @@ class CircuitOutcomeProbabilityArrayLayout(_NicelySerializable):
         The total number of elements in this layout.  In a multi-processor context,
         the number of elements locally owned by the current processor.
 
-    num_elements : int
-        The total number of circuits in this layout.  In a multi-processor context,
-        the number of circuits locally owned by the current processor.
-
     global_layout : CircuitOutcomeProbabilityArrayLayout
         A layout containing all the circuits in their original order, that is the
         same on all processors and doesn't depend on a specific resource allocation.
diff --git a/pygsti/modelmembers/operations/__init__.py b/pygsti/modelmembers/operations/__init__.py
index 579a85523..d0f340c4f 100644
--- a/pygsti/modelmembers/operations/__init__.py
+++ b/pygsti/modelmembers/operations/__init__.py
@@ -42,7 +42,7 @@
 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 """
@@ -86,7 +86,7 @@ def create_from_unitary_mx(unitary_mx, op_type, basis='pp', stdname=None, evotyp
 
                 if op.dim <= 16:  # only do this for up to 2Q operations, otherwise to_dense is too expensive
                     expected_superop_mx = _ot.unitary_to_superop(U, basis)
-                    assert (_np.linalg.norm(op.to_dense('HilbertSchmidt') - expected_superop_mx) < 1e-6), \
+                    assert (_np.linalg.norm(op.to_dense("HilbertSchmidt") - expected_superop_mx) < 1e-6), \
                         "Failure to create Lindblad operation (maybe due the complex log's branch cut?)"
             else:
                 raise ValueError("Unknown operation type '%s'!" % str(typ))
@@ -141,7 +141,7 @@ def create_from_superop_mx(superop_mx, op_type, basis='pp', stdname=None, evotyp
                     ret = ExpErrorgenOp(errorgen)
 
                 if ret.dim <= 16:  # only do this for up to 2Q operations, otherwise to_dense is too expensive
-                    assert(_np.linalg.norm(superop_mx - ret.to_dense('HilbertSchmidt'))
+                    assert(_np.linalg.norm(superop_mx - ret.to_dense("HilbertSchmidt"))
                            < 1e-6), "Failure to create CPTP operation (maybe due the complex log's branch cut?)"
                 return ret
 
@@ -341,11 +341,11 @@ def convert_structure(op):
                 proj_basis = 'PP' if operation.state_space.is_entirely_qubits else basis
                 if ideal_operation == "identity":  # special value
                     postfactor_op = None
-                    error_map_mx = operation.to_dense('HilbertSchmidt')  # error generators are only in HS space
+                    error_map_mx = operation.to_dense("HilbertSchmidt")  # error generators are only in HS space
                 else:
                     postfactor_op = ideal_operation if (ideal_operation is not None) else operation
-                    error_map_mx = _np.dot(operation.to_dense('HilbertSchmidt'),
-                                           _np.linalg.inv(postfactor_op.to_dense('HilbertSchmidt')))
+                    error_map_mx = _np.dot(operation.to_dense("HilbertSchmidt"),
+                                           _np.linalg.inv(postfactor_op.to_dense("HilbertSchmidt")))
 
                 lndtype = LindbladParameterization.cast(to_type)
                 if lndtype.meta == '1+':
@@ -366,7 +366,7 @@ def convert_structure(op):
 
                 if ret.dim <= 81:  # only do this for up to 2-qutrit operations, otherwise to_dense is too expensive
                                    #This should probably be a relative tolerance.
-                    op_diff = _np.linalg.norm(operation.to_dense('HilbertSchmidt') - ret.to_dense('HilbertSchmidt'))
+                    op_diff = _np.linalg.norm(operation.to_dense("HilbertSchmidt") - ret.to_dense("HilbertSchmidt"))
                     assert(op_diff< cptp_truncation_tol), "Failure to create CPTP operation, frobenius norm between original and converted operation large," \
                                             + str(op_diff) + ", (maybe due the complex log's branch cut?)"
                 return ret
diff --git a/pygsti/modelmembers/operations/composederrorgen.py b/pygsti/modelmembers/operations/composederrorgen.py
index f9c321c52..310129993 100644
--- a/pygsti/modelmembers/operations/composederrorgen.py
+++ b/pygsti/modelmembers/operations/composederrorgen.py
@@ -23,7 +23,7 @@
 from pygsti.baseobjs.basis import ExplicitBasis as _ExplicitBasis
 from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel, LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel
 from pygsti.tools import matrixtools as _mt
-
+from pygsti import SpaceT
 
 class ComposedErrorgen(_LinearOperator):
     """
@@ -549,7 +549,7 @@ def remove(self, *factor_indices):
             self.parent._mark_for_rebuild(self)  # of our params may have changed
             self._parent = None  # mark this object for re-allocation
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return this error generator as a sparse matrix
 
@@ -564,7 +564,7 @@ def to_sparse(self, on_space='minimal'):
             mx += eg.to_sparse(on_space)
         return mx
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this error generator as a dense matrix
 
diff --git a/pygsti/modelmembers/operations/composedop.py b/pygsti/modelmembers/operations/composedop.py
index 303b77dc5..d0603cc2c 100644
--- a/pygsti/modelmembers/operations/composedop.py
+++ b/pygsti/modelmembers/operations/composedop.py
@@ -27,6 +27,7 @@
 from pygsti.tools import listtools as _lt
 from pygsti.tools import matrixtools as _mt
 from pygsti.tools import slicetools as _slct
+from pygsti import SpaceT
 
 
 class ComposedOp(_LinearOperator):
@@ -114,9 +115,9 @@ def _update_denserep(self):
         if len(self.factorops) == 0:
             mx = _np.identity(self.state_space.dim, 'd')
         else:
-            mx = self.factorops[0].to_dense(on_space='HilbertSchmidt')
+            mx = self.factorops[0].to_dense("HilbertSchmidt")
             for op in self.factorops[1:]:
-                mx = _np.dot(op.to_dense(on_space='HilbertSchmidt'), mx)
+                mx = _np.dot(op.to_dense("HilbertSchmidt"), mx)
 
         self._rep.base.flags.writeable = True
         self._rep.base[:, :] = mx
@@ -259,7 +260,7 @@ def remove(self, *factorop_indices):
             self.parent._mark_for_rebuild(self)  # of our params may have changed
             self._parent = None  # mark this object for re-allocation
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a sparse matrix
 
@@ -272,7 +273,7 @@ def to_sparse(self, on_space='minimal'):
             mx = op.to_sparse(on_space).dot(mx)
         return mx
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a dense matrix.
 
@@ -395,7 +396,7 @@ def deriv_wrt_params(self, wrt_filter=None):
         numpy array
             Array of derivatives with shape (dimension^2, num_params)
         """
-        typ = complex if any([_np.iscomplexobj(op.to_dense(on_space='minimal'))
+        typ = complex if any([_np.iscomplexobj(op.to_dense("minimal"))
                               for op in self.factorops]) else 'd'
         derivMx = _np.zeros((self.dim, self.dim, self.num_params), typ)
 
@@ -407,16 +408,16 @@ def deriv_wrt_params(self, wrt_filter=None):
             deriv.shape = (self.dim, self.dim, op.num_params)
 
             if i > 0:  # factors before ith
-                pre = self.factorops[0].to_dense(on_space='minimal')
+                pre = self.factorops[0].to_dense("minimal")
                 for opA in self.factorops[1:i]:
-                    pre = _np.dot(opA.to_dense(on_space='minimal'), pre)
+                    pre = _np.dot(opA.to_dense("minimal"), pre)
                 #deriv = _np.einsum("ija,jk->ika", deriv, pre )
                 deriv = _np.transpose(_np.tensordot(deriv, pre, (1, 0)), (0, 2, 1))
 
             if i + 1 < len(self.factorops):  # factors after ith
-                post = self.factorops[i + 1].to_dense(on_space='minimal')
+                post = self.factorops[i + 1].to_dense("minimal")
                 for opA in self.factorops[i + 2:]:
-                    post = _np.dot(opA.to_dense(on_space='minimal'), post)
+                    post = _np.dot(opA.to_dense("minimal"), post)
                 #deriv = _np.einsum("ij,jka->ika", post, deriv )
                 deriv = _np.tensordot(post, deriv, (1, 0))
 
@@ -636,8 +637,8 @@ def transform_inplace(self, s):
             #SPECIAL CASE / HACK: for 1 & 2Q, when holding e^L * T, where T is a static gate
             # then try to gauge transform by setting e^L directly and leaving T alone:
             Smx = s.transform_matrix; Si = s.transform_matrix_inverse
-            Tinv = _np.linalg.inv(self.factorops[0].to_dense(on_space='minimal'))
-            trans_eLT = _np.dot(Si, _np.dot(self.to_dense(on_space='minimal'), Smx))
+            Tinv = _np.linalg.inv(self.factorops[0].to_dense("minimal"))
+            trans_eLT = _np.dot(Si, _np.dot(self.to_dense("minimal"), Smx))
             self.factorops[1].set_dense(_np.dot(trans_eLT, Tinv))  # set_dense(trans_eL)
             return
 
diff --git a/pygsti/modelmembers/operations/denseop.py b/pygsti/modelmembers/operations/denseop.py
index fec40c75e..e5f600182 100644
--- a/pygsti/modelmembers/operations/denseop.py
+++ b/pygsti/modelmembers/operations/denseop.py
@@ -25,7 +25,7 @@
 from pygsti.tools import matrixtools as _mt
 from pygsti.tools import jamiolkowski as _jt
 from pygsti.tools import optools as _ot
-
+from pygsti import SpaceT
 
 def finite_difference_deriv_wrt_params(operation, wrt_filter, eps=1e-7):
     """
@@ -163,7 +163,7 @@ def to_array(self):
         return _np.asarray(self._ptr)
         # *must* be a numpy array for Cython arg conversion
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a sparse matrix.
 
@@ -328,7 +328,7 @@ def _ptr_has_changed(self):
             when the `_ptr` property is changed. """
         self._rep.base_has_changed()
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return the dense array used to represent this operation within its evolution type.
 
@@ -406,7 +406,7 @@ def kraus_operators(self):
         #   -> let reshaped K-th evector be called O_K and dual O_K^dag
         #  = sum_K D_KK O_K rho O_K^dag
         assert(self._basis is not None), "Kraus operator functionality requires specifying a superoperator basis"
-        superop_mx = self.to_dense('HilbertSchmidt'); d = int(_np.round(_np.sqrt(superop_mx.shape[0])))
+        superop_mx = self.to_dense("HilbertSchmidt"); d = int(_np.round(_np.sqrt(superop_mx.shape[0])))
         std_basis = _Basis.cast('std', superop_mx.shape[0])
         choi_mx = _jt.jamiolkowski_iso(superop_mx, self._basis, std_basis) * d  # see NOTE below
         # NOTE: multiply by `d` (density mx dimension) to un-normalize choi_mx as given by
@@ -568,7 +568,7 @@ def _ptr_has_changed(self):
             self._rep.base[:, :] = _ot.unitary_to_superop(self._unitary, self._basis)
         self._rep.base_has_changed()
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return the dense array used to represent this operation within its evolution type.
 
@@ -611,7 +611,7 @@ def to_memoized_dict(self, mmg_memo):
         """
         mm_dict = super().to_memoized_dict(mmg_memo)
 
-        mm_dict['dense_matrix'] = self._encodemx(self.to_dense('Hilbert'))
+        mm_dict['dense_matrix'] = self._encodemx(self.to_dense("Hilbert"))
         mm_dict['basis'] = self._basis.to_nice_serialization()
 
         return mm_dict
@@ -636,7 +636,7 @@ def _is_similar(self, other, rtol, atol):
     @property
     def kraus_operators(self):
         """A list of this operation's Kraus operators as numpy arrays."""
-        return [self.to_dense('Hilbert')]
+        return [self.to_dense("Hilbert")]
 
     def set_kraus_operators(self, kraus_operators):
         """
diff --git a/pygsti/modelmembers/operations/embeddedop.py b/pygsti/modelmembers/operations/embeddedop.py
index 980038fbb..da014dc03 100644
--- a/pygsti/modelmembers/operations/embeddedop.py
+++ b/pygsti/modelmembers/operations/embeddedop.py
@@ -19,8 +19,7 @@
 from pygsti.modelmembers import modelmember as _modelmember
 from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel, LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel
-
-
+from pygsti import SpaceT
 class EmbeddedOp(_LinearOperator):
     """
     An operation containing a single lower (or equal) dimensional operation within it.
@@ -95,7 +94,7 @@ def _create_rep_object(self, evotype, state_space):
     def _update_denserep(self):
         """Performs additional update for the case when we use a dense underlying representation."""
         self._rep.base.flags.writeable = True
-        self._rep.base[:, :] = self.to_dense(on_space='minimal')
+        self._rep.base[:, :] = self.to_dense("minimal")
         self._rep.base.flags.writeable = False
 
     def _update_submember_state_spaces(self, old_parent_state_space, new_parent_state_space):
@@ -240,7 +239,7 @@ def _decomp_op_index(indx, divisors):
                     self._iter_elements_cache[on_space].append(item)
                     yield item
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a sparse matrix
 
@@ -262,7 +261,7 @@ def to_sparse(self, on_space='minimal'):
             finalOp[i, j] = embedded_sparse[gi, gj]
         return finalOp.tocsr()
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a dense matrix
 
diff --git a/pygsti/modelmembers/operations/experrorgenop.py b/pygsti/modelmembers/operations/experrorgenop.py
index aebc9f0c7..9fc420eeb 100644
--- a/pygsti/modelmembers/operations/experrorgenop.py
+++ b/pygsti/modelmembers/operations/experrorgenop.py
@@ -10,7 +10,6 @@
 # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
 #***************************************************************************************************
 
-import warnings as _warnings
 import math
 from typing import Union
 
@@ -25,6 +24,7 @@
 from pygsti.modelmembers import modelmember as _modelmember, term as _term
 from pygsti.modelmembers.errorgencontainer import ErrorGeneratorContainer as _ErrorGeneratorContainer
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
+from pygsti import SpaceT
 
 IMAG_TOL = 1e-7  # tolerance for imaginary part being considered zero
 MAX_EXPONENT = _np.log(_np.finfo('d').max) - 10.0  # so that exp(.) doesn't overflow
@@ -116,7 +116,7 @@ def _update_rep(self, close=False):
         """
         if self._rep_type == 'dense':
             # compute matrix-exponential explicitly
-            self.exp_err_gen = _spl.expm(self.errorgen.to_dense(on_space='HilbertSchmidt'))  # used in deriv_wrt_params
+            self.exp_err_gen = _spl.expm(self.errorgen.to_dense("HilbertSchmidt"))  # used in deriv_wrt_params
 
             dense = self.exp_err_gen
             self._rep.base.flags.writeable = True
@@ -171,7 +171,7 @@ def set_gpindices(self, gpindices, parent, memo=None):
         self.exp_terms_cache = {}
         self.local_term_poly_coeffs = {}
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a dense matrix.
 
@@ -188,7 +188,7 @@ def to_dense(self, on_space='minimal'):
             return _spl.expm(self.errorgen.to_dense(on_space))
 
     #FUTURE: maybe remove this function altogether, as it really shouldn't be called
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a sparse matrix.
 
@@ -240,7 +240,7 @@ def deriv_wrt_params(self, wrt_filter=None):
             #Deriv wrt hamiltonian params
             derrgen = self.errorgen.deriv_wrt_params(None)  # apply filter below; cache *full* deriv
             derrgen.shape = (d2, d2, -1)  # separate 1st d2**2 dim to (d2,d2)
-            dexpL = _d_exp_x(self.errorgen.to_dense(on_space='minimal'), derrgen, self.exp_err_gen)
+            dexpL = _d_exp_x(self.errorgen.to_dense("minimal"), derrgen, self.exp_err_gen)
             derivMx = dexpL.reshape(d2**2, self.num_params)  # [iFlattenedOp,iParam]
 
             assert(_np.linalg.norm(_np.imag(derivMx)) < IMAG_TOL), \
@@ -315,7 +315,7 @@ def hessian_wrt_params(self, wrt_filter1=None, wrt_filter2=None):
             dEdp.shape = (d2, d2, nP)  # separate 1st d2**2 dim to (d2,d2)
             d2Edp2.shape = (d2, d2, nP, nP)  # ditto
 
-            series, series2 = _d2_exp_series(self.errorgen.to_dense(on_space='minimal'), dEdp, d2Edp2)
+            series, series2 = _d2_exp_series(self.errorgen.to_dense("minimal"), dEdp, d2Edp2)
             term1 = series2
             term2 = _np.einsum("ija,jkq->ikaq", series, series)
             d2expL = _np.einsum("ikaq,kj->ijaq", term1 + term2,
@@ -696,7 +696,7 @@ def spam_transform_inplace(self, s, typ):
            or isinstance(s, _gaugegroup.TPSpamGaugeGroupElement):
             U = s.transform_matrix
             Uinv = s.transform_matrix_inverse
-            mx = self.to_dense(on_space='minimal') if self._rep_type == 'dense' else self.to_sparse(on_space='minimal')
+            mx = self.to_dense("minimal") if self._rep_type == 'dense' else self.to_sparse(on_space="minimal")
 
             #just act on postfactor and Lindbladian exponent:
             if typ == "prep":
diff --git a/pygsti/modelmembers/operations/fullcptpop.py b/pygsti/modelmembers/operations/fullcptpop.py
index cf6a1759e..9ea725db3 100644
--- a/pygsti/modelmembers/operations/fullcptpop.py
+++ b/pygsti/modelmembers/operations/fullcptpop.py
@@ -20,6 +20,7 @@
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.tools import jamiolkowski as _jt
 from pygsti.tools import basistools as _bt
+from pygsti import SpaceT
 
 IMAG_TOL = 1e-7
 
@@ -148,7 +149,7 @@ def _update_rep(self):
         else:  # self._reptype == 'dense':
             self._update_dense_rep()
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return the dense array used to represent this operation within its evolution type.
 
diff --git a/pygsti/modelmembers/operations/identitypluserrorgenop.py b/pygsti/modelmembers/operations/identitypluserrorgenop.py
index ec662a9e7..7324d3fef 100644
--- a/pygsti/modelmembers/operations/identitypluserrorgenop.py
+++ b/pygsti/modelmembers/operations/identitypluserrorgenop.py
@@ -23,7 +23,7 @@
 from pygsti.modelmembers.errorgencontainer import ErrorGeneratorContainer as _ErrorGeneratorContainer
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
 from pygsti.tools import matrixtools as _mt
-
+from pygsti import SpaceT
 IMAG_TOL = 1e-7  # tolerance for imaginary part being considered zero
 TODENSE_TRUNCATE = 1e-11
 
@@ -106,7 +106,7 @@ def _update_rep(self, close=False):
         if self._rep_type == 'dense':
             # compute 1 + errorgen explicitly
             self._rep.base.flags.writeable = True
-            self._rep.base[:, :] = self._ident + self.errorgen.to_dense(on_space='HilbertSchmidt')
+            self._rep.base[:, :] = self._ident + self.errorgen.to_dense("HilbertSchmidt")
             self._rep.base.flags.writeable = False
         else:
             pass  # nothing to do -- no need to even notify OpRepIdentityPlusErrorgen that errorgen has changed
@@ -134,7 +134,7 @@ def set_gpindices(self, gpindices, parent, memo=None):
         self.terms = {}  # clear terms cache since param indices have changed now
         self.local_term_poly_coeffs = {}
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a dense matrix.
 
@@ -152,7 +152,7 @@ def to_dense(self, on_space='minimal'):
             # Note: above fails with shape mismatch if try to get dense Hilbert-space op - is this desired?
 
     #FUTURE: maybe remove this function altogether, as it really shouldn't be called
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a sparse matrix.
 
diff --git a/pygsti/modelmembers/operations/lindbladerrorgen.py b/pygsti/modelmembers/operations/lindbladerrorgen.py
index cbefec4c8..40b1b478e 100644
--- a/pygsti/modelmembers/operations/lindbladerrorgen.py
+++ b/pygsti/modelmembers/operations/lindbladerrorgen.py
@@ -30,6 +30,7 @@
 from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel
 from pygsti.tools import matrixtools as _mt
 from pygsti.tools import optools as _ot
+from pygsti import SpaceT
 
 IMAG_TOL = 1e-7  # tolerance for imaginary part being considered zero
 
@@ -701,7 +702,7 @@ def _update_rep(self):
 
         self._onenorm_upbound = onenorm
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this error generator as a dense matrix.
 
@@ -731,7 +732,7 @@ def to_dense(self, on_space='minimal'):
         else:  # dense rep
             return self._rep.to_dense(on_space)
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the error generator as a sparse matrix.
 
diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py
index 7cd00c073..a2ab0c0ea 100644
--- a/pygsti/modelmembers/operations/linearop.py
+++ b/pygsti/modelmembers/operations/linearop.py
@@ -15,7 +15,7 @@
 from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex
 from pygsti.modelmembers import modelmember as _modelmember
 from pygsti.tools import optools as _ot
-
+from pygsti import SpaceT
 
 #Note on initialization sequence of Operations within a Model:
 # 1) a Model is constructed (empty)
@@ -128,7 +128,7 @@ def set_time(self, t):
         """
         pass
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a dense matrix.
 
@@ -170,7 +170,7 @@ def acton(self, state, on_space='minimal'):
         #Build a State around output_rep
         return _state.StaticState(output_rep.to_dense(on_space), None, self._evotype, self.state_space)
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a sparse matrix.
 
@@ -406,11 +406,11 @@ def frobeniusdist_squared(self, other_op, transform=None, inv_transform=None):
         float
         """
         if transform is None and inv_transform is None:
-            return _ot.frobeniusdist_squared(self.to_dense(on_space='minimal'), other_op.to_dense(on_space='minimal'))
+            return _ot.frobeniusdist_squared(self.to_dense("minimal"), other_op.to_dense("minimal"))
         else:
             return _ot.frobeniusdist_squared(_np.dot(
-                inv_transform, _np.dot(self.to_dense(on_space='minimal'), transform)),
-                other_op.to_dense(on_space='minimal'))
+                inv_transform, _np.dot(self.to_dense("minimal"), transform)),
+                other_op.to_dense("minimal"))
 
     def frobeniusdist(self, other_op, transform=None, inv_transform=None):
         """
@@ -461,13 +461,13 @@ def residuals(self, other_op, transform=None, inv_transform=None):
         numpy.ndarray
             A 1D-array of size equal to that of the flattened operation matrix.
         """
-        dense_self = self.to_dense(on_space='minimal')
+        dense_self = self.to_dense("minimal")
         if transform is not None:
             assert inv_transform is not None
             dense_self = inv_transform @ (dense_self @ transform)
         else:
             assert inv_transform is None
-        return (dense_self - other_op.to_dense(on_space='minimal')).ravel()
+        return (dense_self - other_op.to_dense("minimal")).ravel()
 
 
     def jtracedist(self, other_op, transform=None, inv_transform=None):
@@ -493,11 +493,11 @@ def jtracedist(self, other_op, transform=None, inv_transform=None):
         float
         """
         if transform is None and inv_transform is None:
-            return _ot.jtracedist(self.to_dense(on_space='minimal'), other_op.to_dense(on_space='minimal'))
+            return _ot.jtracedist(self.to_dense("minimal"), other_op.to_dense("minimal"))
         else:
             return _ot.jtracedist(_np.dot(
-                inv_transform, _np.dot(self.to_dense(on_space='minimal'), transform)),
-                other_op.to_dense(on_space='minimal'))
+                inv_transform, _np.dot(self.to_dense("minimal"), transform)),
+                other_op.to_dense("minimal"))
 
     def diamonddist(self, other_op, transform=None, inv_transform=None):
         """
@@ -522,11 +522,11 @@ def diamonddist(self, other_op, transform=None, inv_transform=None):
         float
         """
         if transform is None and inv_transform is None:
-            return _ot.diamonddist(self.to_dense(on_space='minimal'), other_op.to_dense(on_space='minimal'))
+            return _ot.diamonddist(self.to_dense("minimal"), other_op.to_dense("minimal"))
         else:
             return _ot.diamonddist(_np.dot(
-                inv_transform, _np.dot(self.to_dense(on_space='minimal'), transform)),
-                other_op.to_dense(on_space='minimal'))
+                inv_transform, _np.dot(self.to_dense("minimal"), transform)),
+                other_op.to_dense("minimal"))
 
     def transform_inplace(self, s):
         """
@@ -552,7 +552,7 @@ def transform_inplace(self, s):
         """
         Smx = s.transform_matrix
         Si = s.transform_matrix_inverse
-        self.set_dense(_np.dot(Si, _np.dot(self.to_dense(on_space='minimal'), Smx)))
+        self.set_dense(_np.dot(Si, _np.dot(self.to_dense("minimal"), Smx)))
 
     def spam_transform_inplace(self, s, typ):
         """
@@ -580,9 +580,9 @@ def spam_transform_inplace(self, s, typ):
         None
         """
         if typ == 'prep':
-            self.set_dense(_np.dot(s.transform_matrix_inverse, self.to_dense(on_space='minimal')))
+            self.set_dense(_np.dot(s.transform_matrix_inverse, self.to_dense("minimal")))
         elif typ == 'effect':
-            self.set_dense(_np.dot(self.to_dense(on_space='minimal'), s.transform_matrix))
+            self.set_dense(_np.dot(self.to_dense("minimal"), s.transform_matrix))
         else:
             raise ValueError("Invalid `typ` argument: %s" % typ)
 
@@ -616,7 +616,7 @@ def depolarize(self, amount):
         else:
             assert(len(amount) == self.dim - 1)
             D = _np.diag([1] + list(1.0 - _np.array(amount, 'd')))
-        self.set_dense(_np.dot(D, self.to_dense(on_space='minimal')))
+        self.set_dense(_np.dot(D, self.to_dense("minimal")))
 
     def rotate(self, amount, mx_basis="gm"):
         """
@@ -647,7 +647,7 @@ def rotate(self, amount, mx_basis="gm"):
         None
         """
         rotnMx = _ot.rotation_gate_mx(amount, mx_basis)
-        self.set_dense(_np.dot(rotnMx, self.to_dense(on_space='minimal')))
+        self.set_dense(_np.dot(rotnMx, self.to_dense("minimal")))
 
     def deriv_wrt_params(self, wrt_filter=None):
         """
@@ -812,7 +812,7 @@ def finite_difference_deriv_wrt_params(operation, wrt_filter, eps=1e-7):
         An M by N matrix where M is the number of operation elements and
         N is the number of operation parameters.
     """
-    dense_operation = operation.to_dense(on_space='minimal')
+    dense_operation = operation.to_dense("minimal")
     dim = dense_operation.shape[0]
     #operation.from_vector(operation.to_vector()) #ensure we call from_vector w/close=False first
     op2 = operation.copy()
@@ -826,7 +826,7 @@ def finite_difference_deriv_wrt_params(operation, wrt_filter, eps=1e-7):
         p_plus_dp = p.copy()
         p_plus_dp[i] += eps
         op2.from_vector(p_plus_dp)
-        fd_deriv[:, :, ii] = (op2.to_dense(on_space='minimal') - dense_operation) / eps
+        fd_deriv[:, :, ii] = (op2.to_dense("minimal") - dense_operation) / eps
 
     fd_deriv.shape = [dim**2, len(wrt_filter)]
     return fd_deriv
@@ -865,7 +865,7 @@ def finite_difference_hessian_wrt_params(operation, wrt_filter1, wrt_filter2, ep
         N1 and N2 are numbers of operation parameters.
     """
     #operation.from_vector(operation.to_vector()) #ensure we call from_vector w/close=False first
-    dense_operation = operation.to_dense(on_space='minimal')
+    dense_operation = operation.to_dense("minimal")
     fd_deriv0 = finite_difference_deriv_wrt_params(operation, wrt_filter1, eps=eps)
 
     if wrt_filter2 is None:
diff --git a/pygsti/modelmembers/operations/repeatedop.py b/pygsti/modelmembers/operations/repeatedop.py
index e9a0bdca3..60cccc361 100644
--- a/pygsti/modelmembers/operations/repeatedop.py
+++ b/pygsti/modelmembers/operations/repeatedop.py
@@ -15,6 +15,7 @@
 
 from pygsti.modelmembers.operations.linearop import LinearOperator as _LinearOperator
 from pygsti.evotypes import Evotype as _Evotype
+from pygsti import SpaceT
 
 
 class RepeatedOp(_LinearOperator):
@@ -77,7 +78,7 @@ def set_time(self, t):
         """
         self.repeated_op.set_time(t)
 
-    def to_sparse(self, on_space='minimal'):
+    def to_sparse(self, on_space: SpaceT='minimal'):
         """
         Return the operation as a sparse matrix
 
@@ -94,7 +95,7 @@ def to_sparse(self, on_space='minimal'):
             mx = mx.dot(op)
         return mx
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a dense matrix.
 
@@ -192,7 +193,7 @@ def deriv_wrt_params(self, wrt_filter=None):
         numpy array
             Array of derivatives with shape (dimension^2, num_params)
         """
-        mx = self.repeated_op.to_dense(on_space='minimal')
+        mx = self.repeated_op.to_dense("minimal")
 
         mx_powers = {0: _np.identity(self.dim, 'd'), 1: mx}
         for i in range(2, self.num_repetitions):
diff --git a/pygsti/modelmembers/operations/staticcliffordop.py b/pygsti/modelmembers/operations/staticcliffordop.py
index eefd24918..33ff6b850 100644
--- a/pygsti/modelmembers/operations/staticcliffordop.py
+++ b/pygsti/modelmembers/operations/staticcliffordop.py
@@ -19,7 +19,7 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
-
+from pygsti import SpaceT
 
 class StaticCliffordOp(_LinearOperator, _NoErrorGeneratorInterface):
     """
@@ -168,7 +168,7 @@ def total_term_magnitude_deriv(self):
         """
         return _np.empty((0,), 'd')
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return the dense array used to represent this operation within its evolution type.
 
diff --git a/pygsti/modelmembers/operations/staticstdop.py b/pygsti/modelmembers/operations/staticstdop.py
index 77065ab89..3bf63a405 100644
--- a/pygsti/modelmembers/operations/staticstdop.py
+++ b/pygsti/modelmembers/operations/staticstdop.py
@@ -19,6 +19,7 @@
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
 from pygsti.tools import internalgates as _itgs
+from pygsti import SpaceT
 
 
 class StaticStandardOp(_LinearOperator, _NoErrorGeneratorInterface):
@@ -58,7 +59,7 @@ def __init__(self, name, basis='pp', evotype="default", state_space=None):
         rep = evotype.create_standard_rep(name, basis, state_space)
         _LinearOperator.__init__(self, rep, evotype)
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return the dense array used to represent this operation within its evolution type.
 
diff --git a/pygsti/modelmembers/operations/stochasticop.py b/pygsti/modelmembers/operations/stochasticop.py
index f863f05ac..8c331e3ad 100644
--- a/pygsti/modelmembers/operations/stochasticop.py
+++ b/pygsti/modelmembers/operations/stochasticop.py
@@ -19,7 +19,7 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
-
+from pygsti import SpaceT
 
 class StochasticNoiseOp(_LinearOperator, _KrausOperatorInterface):
     """
@@ -100,7 +100,7 @@ def _get_rate_poly_dicts(self):
             keys of dicts <=> poly terms, e.g. (1,1) <=> x1^2) """
         return [{(i, i): 1.0} for i in range(self.basis.size - 1)]  # rates are just parameters squared
 
-    def to_dense(self, on_space='minimal'):
+    def to_dense(self, on_space: SpaceT='minimal'):
         """
         Return this operation as a dense matrix.
 
@@ -290,7 +290,7 @@ def total_term_magnitude_deriv(self):
     @property
     def kraus_operators(self):
         """A list of this operation's Kraus operators as numpy arrays."""
-        kraus_ops = [_np.sqrt(urate) * urep.to_dense('Hilbert')
+        kraus_ops = [_np.sqrt(urate) * urep.to_dense("Hilbert")
                      for urate, urep in zip(self._rep.unitary_rates, self._rep.unitary_reps)]
         return kraus_ops
 
diff --git a/pygsti/modelmembers/povms/composedeffect.py b/pygsti/modelmembers/povms/composedeffect.py
index bd81be38e..c43d08d0c 100644
--- a/pygsti/modelmembers/povms/composedeffect.py
+++ b/pygsti/modelmembers/povms/composedeffect.py
@@ -16,7 +16,7 @@
 from pygsti.modelmembers.povms.effect import POVMEffect as _POVMEffect
 from pygsti.modelmembers import modelmember as _modelmember, term as _term
 from pygsti.modelmembers.states.staticstate import StaticState as _StaticState
-
+from pygsti import SpaceT
 
 class ComposedPOVMEffect(_POVMEffect):  # , _ErrorMapContainer
     """
@@ -113,7 +113,7 @@ def set_gpindices(self, gpindices, parent, memo=None):
         self.local_term_poly_coeffs = {}
         super().set_gpindices(gpindices, parent, memo)
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this POVM effect vector as a (dense) numpy array.
 
@@ -299,7 +299,7 @@ def deriv_wrt_params(self, wrt_filter=None):
         numpy array
             Array of derivatives, shape == (dimension, num_params)
         """
-        dmVec = self.effect_vec.to_dense(on_space='minimal')
+        dmVec = self.effect_vec.to_dense("minimal")
 
         derrgen = self.error_map.deriv_wrt_params(wrt_filter)  # shape (dim*dim, n_params)
         derrgen.shape = (self.dim, self.dim, derrgen.shape[1])  # => (dim,dim,n_params)
@@ -332,7 +332,7 @@ def hessian_wrt_params(self, wrt_filter1=None, wrt_filter2=None):
         numpy array
             Hessian with shape (dimension, num_params1, num_params2)
         """
-        dmVec = self.effect_vec.to_dense(on_space='minimal')
+        dmVec = self.effect_vec.to_dense("minimal")
 
         herrgen = self.error_map.hessian_wrt_params(wrt_filter1, wrt_filter2)  # shape (dim*dim, nParams1, nParams2)
         herrgen.shape = (self.dim, self.dim, herrgen.shape[1], herrgen.shape[2])  # => (dim,dim,nParams1, nParams2)
diff --git a/pygsti/modelmembers/povms/computationaleffect.py b/pygsti/modelmembers/povms/computationaleffect.py
index 96502aeac..151170ef2 100644
--- a/pygsti/modelmembers/povms/computationaleffect.py
+++ b/pygsti/modelmembers/povms/computationaleffect.py
@@ -22,6 +22,7 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
+from pygsti import SpaceT
 
 try:
     from pygsti.tools import fastcalc as _fastcalc
@@ -191,7 +192,7 @@ def _is_similar(self, other, rtol, atol):
         return (_np.array_equal(self._rep.zvals, other._rep.zvals)
                 and self._rep.basis == other._rep.basis)
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this POVM effect vector as a (dense) numpy array.
 
diff --git a/pygsti/modelmembers/povms/conjugatedeffect.py b/pygsti/modelmembers/povms/conjugatedeffect.py
index df2745e33..dcfaad872 100644
--- a/pygsti/modelmembers/povms/conjugatedeffect.py
+++ b/pygsti/modelmembers/povms/conjugatedeffect.py
@@ -16,6 +16,7 @@
 from pygsti.modelmembers.povms.effect import POVMEffect as _POVMEffect
 from pygsti.modelmembers import term as _term
 from pygsti.tools import matrixtools as _mt
+from pygsti import SpaceT
 
 
 class DenseEffectInterface(object):
@@ -161,7 +162,7 @@ def parameter_labels(self):
         """
         return self.state.parameter_labels
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this POVM effect vector as a (dense) numpy array.
 
@@ -200,7 +201,7 @@ def hilbert_schmidt_size(self):
 
     def __str__(self):
         s = "%s with dimension %d\n" % (self.__class__.__name__, self.dim)
-        s += _mt.mx_to_string(self.to_dense(on_space='minimal'), width=4, prec=2)
+        s += _mt.mx_to_string(self.to_dense("minimal"), width=4, prec=2)
         return s
 
     def submembers(self):
diff --git a/pygsti/modelmembers/povms/tensorprodeffect.py b/pygsti/modelmembers/povms/tensorprodeffect.py
index b1a68b5ad..0c2acf4d9 100644
--- a/pygsti/modelmembers/povms/tensorprodeffect.py
+++ b/pygsti/modelmembers/povms/tensorprodeffect.py
@@ -21,6 +21,7 @@
 from pygsti.tools import listtools as _lt
 from pygsti.tools import matrixtools as _mt
 from pygsti.tools import slicetools as _slct
+from pygsti import SpaceT
 
 
 class TensorProductPOVMEffect(_POVMEffect):
@@ -113,7 +114,7 @@ def parameter_labels(self):
             vl[povm_local_inds] = povm.parameter_labels
         return vl
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this POVM effect vector as a (dense) numpy array.
 
@@ -206,13 +207,13 @@ def taylor_order_terms(self, order, max_polynomial_vars=100, return_coeff_polys=
                 # as they'd only be used if the evotype was like densitymx and needed to convert to
                 # a dense superoperator.
                 pre_state_rep = self._evotype.create_tensorproduct_state_rep(
-                    [self._evotype.create_pure_state_rep(f.pre_effect.to_dense(on_space='Hilbert'), None,
+                    [self._evotype.create_pure_state_rep(f.pre_effect.to_dense("Hilbert"), None,
                                                          f.pre_effect.state_space)
                      for f in factors if (f.pre_effect is not None)], self.state_space)
                 pre_rep = self._evotype.create_conjugatedstate_effect_rep(pre_state_rep)
 
                 post_state_rep = self._evotype.create_tensorproduct_state_rep(
-                    [self._evotype.create_pure_state_rep(f.post_effect.to_dense(on_space='Hilbert'), None,
+                    [self._evotype.create_pure_state_rep(f.post_effect.to_dense("Hilbert"), None,
                                                          f.post_effect.state_space)
                      for f in factors if (f.post_effect is not None)], self.state_space)
                 post_rep = self._evotype.create_conjugatedstate_effect_rep(post_state_rep)
@@ -332,11 +333,11 @@ def deriv_wrt_params(self, wrt_filter=None):
         numpy array
             Array of derivatives, shape == (dimension, num_params)
         """
-        typ = self.factors[0].to_dense(on_space='minimal').dtype if len(self.factors) > 0 else 'd'
+        typ = self.factors[0].to_dense("minimal").dtype if len(self.factors) > 0 else 'd'
 
         #HACK to deal with fact that output of to_dense is really what is differentiated
         # but this may not match self.dim == self.state_space.dim, e.g. for pure state vecs.
-        dims = [len(fct.to_dense(on_space='minimal')) for fct in self.factors]
+        dims = [len(fct.to_dense("minimal")) for fct in self.factors]
         dim = int(_np.prod(dims))
 
         derivMx = _np.zeros((dim, self.num_params), typ)
@@ -351,15 +352,15 @@ def deriv_wrt_params(self, wrt_filter=None):
             deriv.shape = (fct_dim, vec.num_params)
 
             if i > 0:  # factors before ith
-                pre = self.factors[0][self.effectLbls[0]].to_dense(on_space='minimal')
+                pre = self.factors[0][self.effectLbls[0]].to_dense("minimal")
                 for j, fctA in enumerate(self.factors[1:i], start=1):
-                    pre = _np.kron(pre, fctA[self.effectLbls[j]].to_dense(on_space='minimal'))
+                    pre = _np.kron(pre, fctA[self.effectLbls[j]].to_dense("minimal"))
                 deriv = _np.kron(pre[:, None], deriv)  # add a dummy 1-dim to 'pre' and do kron properly...
 
             if i + 1 < len(self.factors):  # factors after ith
-                post = self.factors[i + 1][self.effectLbls[i + 1]].to_dense(on_space='minimal')
+                post = self.factors[i + 1][self.effectLbls[i + 1]].to_dense("minimal")
                 for j, fctA in enumerate(self.factors[i + 2:], start=i + 2):
-                    post = _np.kron(post, fctA[self.effectLbls[j]].to_dense(on_space='minimal'))
+                    post = _np.kron(post, fctA[self.effectLbls[j]].to_dense("minimal"))
                 deriv = _np.kron(deriv, post[:, None])  # add a dummy 1-dim to 'post' and do kron properly...
 
             assert(fct_local_inds is not None), \
@@ -388,6 +389,6 @@ def __str__(self):
         #s += _mt.mx_to_string(ar, width=4, prec=2)
 
         # factors are POVMs
-        s += " x ".join([_mt.mx_to_string(fct[self.effectLbls[i]].to_dense(on_space='minimal'), width=4, prec=2)
+        s += " x ".join([_mt.mx_to_string(fct[self.effectLbls[i]].to_dense("minimal"), width=4, prec=2)
                          for i, fct in enumerate(self.factors)])
         return s
diff --git a/pygsti/modelmembers/states/composedstate.py b/pygsti/modelmembers/states/composedstate.py
index 8b62c8aef..fd5af4220 100644
--- a/pygsti/modelmembers/states/composedstate.py
+++ b/pygsti/modelmembers/states/composedstate.py
@@ -17,7 +17,7 @@
 from pygsti.modelmembers.states.staticstate import StaticState as _StaticState
 from pygsti.modelmembers import modelmember as _modelmember, term as _term
 from pygsti.modelmembers.errorgencontainer import ErrorMapContainer as _ErrorMapContainer
-
+from pygsti import SpaceT
 
 class ComposedState(_State):  # , _ErrorMapContainer
     """
@@ -115,7 +115,7 @@ def set_gpindices(self, gpindices, parent, memo=None):
         self.local_term_poly_coeffs = {}
         super().set_gpindices(gpindices, parent, memo)
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this state vector as a (dense) numpy array.
 
@@ -291,7 +291,7 @@ def deriv_wrt_params(self, wrt_filter=None):
         numpy array
             Array of derivatives, shape == (dimension, num_params)
         """
-        dmVec = self.state_vec.to_dense(on_space='minimal')
+        dmVec = self.state_vec.to_dense("minimal")
 
         derrgen = self.error_map.deriv_wrt_params(wrt_filter)  # shape (dim*dim, n_params)
         derrgen.shape = (self.dim, self.dim, derrgen.shape[1])  # => (dim,dim,n_params)
@@ -323,7 +323,7 @@ def hessian_wrt_params(self, wrt_filter1=None, wrt_filter2=None):
         numpy array
             Hessian with shape (dimension, num_params1, num_params2)
         """
-        dmVec = self.state_vec.to_dense(on_space='minimal')
+        dmVec = self.state_vec.to_dense("minimal")
 
         herrgen = self.error_map.hessian_wrt_params(wrt_filter1, wrt_filter2)  # shape (dim*dim, nParams1, nParams2)
         herrgen.shape = (self.dim, self.dim, herrgen.shape[1], herrgen.shape[2])  # => (dim,dim,nParams1, nParams2)
diff --git a/pygsti/modelmembers/states/computationalstate.py b/pygsti/modelmembers/states/computationalstate.py
index dc23666e0..a27cc07fb 100644
--- a/pygsti/modelmembers/states/computationalstate.py
+++ b/pygsti/modelmembers/states/computationalstate.py
@@ -22,6 +22,7 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
+from pygsti import SpaceT
 
 try:
     from pygsti.tools import fastcalc as _fastcalc
@@ -145,7 +146,7 @@ def __init__(self, zvals, basis='pp', evotype="default", state_space=None):
         rep = evotype.create_computational_state_rep(self._zvals, basis, state_space)
         _State.__init__(self, rep, evotype)
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this state vector as a (dense) numpy array.
 
@@ -167,8 +168,8 @@ def to_dense(self, on_space='minimal', scratch=None):
         numpy.ndarray
         """
         from .staticpurestate import StaticPureState as _StaticPureState
-        v0 = _StaticPureState(_np.array((1, 0), complex), basis='pp', evotype=self._evotype).to_dense('minimal')
-        v1 = _StaticPureState(_np.array((0, 1), complex), basis='pp', evotype=self._evotype).to_dense('minimal')
+        v0 = _StaticPureState(_np.array((1, 0), complex), basis='pp', evotype=self._evotype).to_dense("minimal")
+        v1 = _StaticPureState(_np.array((0, 1), complex), basis='pp', evotype=self._evotype).to_dense("minimal")
         factor_dim = len(v0)
         v = (v0, v1)
 
diff --git a/pygsti/modelmembers/states/densestate.py b/pygsti/modelmembers/states/densestate.py
index b645f7795..714b39913 100644
--- a/pygsti/modelmembers/states/densestate.py
+++ b/pygsti/modelmembers/states/densestate.py
@@ -21,6 +21,7 @@
 from pygsti.tools import basistools as _bt
 from pygsti.tools import matrixtools as _mt
 from pygsti.tools import optools as _ot
+from pygsti import SpaceT
 
 
 class DenseStateInterface(object):
@@ -134,7 +135,7 @@ def __complex__(self): return complex(self.columnvec)
 
     def __str__(self):
         s = "%s with dimension %d\n" % (self.__class__.__name__, self.dim)
-        s += _mt.mx_to_string(self.to_dense(on_space='minimal'), width=4, prec=2)
+        s += _mt.mx_to_string(self.to_dense("minimal"), width=4, prec=2)
         return s
 
 
@@ -189,7 +190,7 @@ def _ptr_has_changed(self):
         """
         self._rep.base_has_changed()
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return the dense array used to represent this state within its evolution type.
 
@@ -304,7 +305,7 @@ def _ptr_has_changed(self):
             self._rep.base[:] = _bt.change_basis(_ot.state_to_dmvec(self._purevec), 'std', self._basis)
         self._rep.base_has_changed()
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return the dense array used to represent this state within its evolution type.
 
@@ -351,7 +352,7 @@ def to_memoized_dict(self, mmg_memo):
         """
         mm_dict = super().to_memoized_dict(mmg_memo)
 
-        mm_dict['dense_state_vector'] = self._encodemx(self.to_dense('Hilbert'))
+        mm_dict['dense_state_vector'] = self._encodemx(self.to_dense("Hilbert"))
         mm_dict['basis'] = self._basis.to_nice_serialization() if (self._basis is not None) else None
         return mm_dict
 
diff --git a/pygsti/modelmembers/states/purestate.py b/pygsti/modelmembers/states/purestate.py
index a319a0a70..6986bfe6d 100644
--- a/pygsti/modelmembers/states/purestate.py
+++ b/pygsti/modelmembers/states/purestate.py
@@ -20,6 +20,7 @@
 from pygsti.baseobjs.polynomial import Polynomial as _Polynomial
 from pygsti.tools import basistools as _bt
 from pygsti.tools import optools as _ot
+from pygsti import SpaceT
 
 
 #TODO: figure out what to do with this class when we wire up term calcs??
@@ -86,7 +87,7 @@ def __init__(self, pure_state, evotype='default', dm_basis='pp'):
         #_State.__init__(self, rep, evotype)
         #self.init_gpindices()  # initialize our gpindices based on sub-members
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this state vector as a (dense) numpy array.
 
@@ -108,7 +109,7 @@ def to_dense(self, on_space='minimal', scratch=None):
         numpy.ndarray
         """
         assert(on_space in ('minimal', 'HilbertSchmidt'))
-        dmVec_std = _ot.state_to_dmvec(self.pure_state.to_dense(on_space='Hilbert'))
+        dmVec_std = _ot.state_to_dmvec(self.pure_state.to_dense("Hilbert"))
         return _bt.change_basis(dmVec_std, 'std', self.basis)
 
     def taylor_order_terms(self, order, max_polynomial_vars=100, return_coeff_polys=False):
diff --git a/pygsti/modelmembers/states/state.py b/pygsti/modelmembers/states/state.py
index 70024d5ce..da693a957 100644
--- a/pygsti/modelmembers/states/state.py
+++ b/pygsti/modelmembers/states/state.py
@@ -17,6 +17,7 @@
 from pygsti.modelmembers import modelmember as _modelmember
 from pygsti.baseobjs import _compatibility as _compat
 from pygsti.tools import optools as _ot
+from pygsti import SpaceT
 
 
 class State(_modelmember.ModelMember):
@@ -104,7 +105,7 @@ def set_time(self, t):
         """
         pass
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this state vector as a (dense) numpy array.
 
@@ -320,12 +321,12 @@ def frobeniusdist_squared(self, other_spam_vec, transform=None,
         -------
         float
         """
-        vec = self.to_dense(on_space='minimal')
+        vec = self.to_dense("minimal")
         if inv_transform is None:
-            return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense(on_space='minimal'))
+            return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense("minimal"))
         else:
             return _ot.frobeniusdist_squared(_np.dot(inv_transform, vec),
-                                             other_spam_vec.to_dense(on_space='minimal'))
+                                             other_spam_vec.to_dense("minimal"))
 
     def residuals(self, other_spam_vec, transform=None, inv_transform=None):
         """
@@ -349,10 +350,10 @@ def residuals(self, other_spam_vec, transform=None, inv_transform=None):
         -------
         float
         """
-        vec = self.to_dense(on_space='minimal')
+        vec = self.to_dense("minimal")
         if inv_transform is not None:
             vec = inv_transform @ vec
-        return (vec - other_spam_vec.to_dense(on_space='minimal')).ravel()
+        return (vec - other_spam_vec.to_dense("minimal")).ravel()
 
 
     def transform_inplace(self, s):
@@ -378,7 +379,7 @@ def transform_inplace(self, s):
         None
         """
         Si = s.transform_matrix_inverse
-        self.set_dense(_np.dot(Si, self.to_dense(on_space='minimal')))
+        self.set_dense(_np.dot(Si, self.to_dense("minimal")))
 
     def depolarize(self, amount):
         """
@@ -407,7 +408,7 @@ def depolarize(self, amount):
         else:
             assert(len(amount) == self.dim - 1)
             D = _np.diag([1] + list(1.0 - _np.array(amount, 'd')))
-        self.set_dense(_np.dot(D, self.to_dense(on_space='minimal')))
+        self.set_dense(_np.dot(D, self.to_dense("minimal")))
 
     @property
     def num_params(self):
@@ -541,7 +542,7 @@ def _to_vector(v):
         numpy array
         """
         if isinstance(v, State):
-            vector = v.to_dense(on_space='minimal').copy()
+            vector = v.to_dense("minimal").copy()
             vector.shape = (vector.size, 1)
         elif isinstance(v, _np.ndarray):
             vector = v.copy()
diff --git a/pygsti/modelmembers/states/tensorprodstate.py b/pygsti/modelmembers/states/tensorprodstate.py
index 894b76df8..bb4550cb4 100644
--- a/pygsti/modelmembers/states/tensorprodstate.py
+++ b/pygsti/modelmembers/states/tensorprodstate.py
@@ -21,7 +21,7 @@
 from pygsti.baseobjs import statespace as _statespace
 from pygsti.tools import listtools as _lt
 from pygsti.tools import matrixtools as _mt
-
+from pygsti import SpaceT
 
 class TensorProductState(_State):
     """
@@ -79,7 +79,7 @@ def parameter_labels(self):
             vl[factor_local_inds] = factor_state.parameter_labels
         return vl
 
-    def to_dense(self, on_space='minimal', scratch=None):
+    def to_dense(self, on_space: SpaceT='minimal', scratch=None):
         """
         Return this state vector as a (dense) numpy array.
 
@@ -280,11 +280,11 @@ def deriv_wrt_params(self, wrt_filter=None):
         numpy array
             Array of derivatives, shape == (dimension, num_params)
         """
-        typ = self.factors[0].to_dense(on_space='minimal').dtype if len(self.factors) > 0 else 'd'
+        typ = self.factors[0].to_dense("minimal").dtype if len(self.factors) > 0 else 'd'
 
         #HACK to deal with fact that output of to_dense is really what is differentiated
         # but this may not match self.dim == self.state_space.dim, e.g. for pure state vecs.
-        dims = [len(fct.to_dense(on_space='minimal')) for fct in self.factors]
+        dims = [len(fct.to_dense("minimal")) for fct in self.factors]
         dim = int(_np.prod(dims))
 
         derivMx = _np.zeros((dim, self.num_params), typ)
@@ -299,15 +299,15 @@ def deriv_wrt_params(self, wrt_filter=None):
             deriv.shape = (fct_dim, vec.num_params)
 
             if i > 0:  # factors before ith
-                pre = self.factors[0].to_dense(on_space='minimal')
+                pre = self.factors[0].to_dense("minimal")
                 for vecA in self.factors[1:i]:
-                    pre = _np.kron(pre, vecA.to_dense(on_space='minimal'))
+                    pre = _np.kron(pre, vecA.to_dense("minimal"))
                 deriv = _np.kron(pre[:, None], deriv)  # add a dummy 1-dim to 'pre' and do kron properly...
 
             if i + 1 < len(self.factors):  # factors after ith
-                post = self.factors[i + 1].to_dense(on_space='minimal')
+                post = self.factors[i + 1].to_dense("minimal")
                 for vecA in self.factors[i + 2:]:
-                    post = _np.kron(post, vecA.to_dense(on_space='minimal'))
+                    post = _np.kron(post, vecA.to_dense("minimal"))
                 deriv = _np.kron(deriv, post[:, None])  # add a dummy 1-dim to 'post' and do kron properly...
 
             assert(fct_local_inds is not None), \
@@ -336,5 +336,5 @@ def __str__(self):
         #s += _mt.mx_to_string(ar, width=4, prec=2)
 
         # factors are just other States
-        s += " x ".join([_mt.mx_to_string(fct.to_dense(on_space='minimal'), width=4, prec=2) for fct in self.factors])
+        s += " x ".join([_mt.mx_to_string(fct.to_dense("minimal"), width=4, prec=2) for fct in self.factors])
         return s
diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py
index 01eed989b..0fc571bdf 100644
--- a/pygsti/models/explicitmodel.py
+++ b/pygsti/models/explicitmodel.py
@@ -43,6 +43,7 @@
 from pygsti.tools import fogitools as _fogit
 from pygsti.tools import slicetools as _slct
 from pygsti.tools import listtools as _lt
+from pygsti import SpaceT
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
 
 
@@ -1529,7 +1530,7 @@ def add_availability(opkey, op):
                     #    observed_sslbls.update(sslbls)
 
                 if gn not in gate_unitaries or gate_unitaries[gn] is None:
-                    U = _ot.superop_to_unitary(op.to_dense('HilbertSchmidt'), self.basis) \
+                    U = _ot.superop_to_unitary(op.to_dense("HilbertSchmidt"), self.basis) \
                         if (op is not None) else None  # U == None indicates "unknown, up until this point"
 
                     Ulocal = extract_unitary(U, all_sslbls, sslbls)
diff --git a/pygsti/models/gaugegroup.py b/pygsti/models/gaugegroup.py
index 9f8c83a1f..d42390808 100644
--- a/pygsti/models/gaugegroup.py
+++ b/pygsti/models/gaugegroup.py
@@ -20,7 +20,7 @@
 from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 from pygsti.evotypes.evotype import Evotype as _Evotype
 from pygsti.tools.optools import superop_to_unitary, unitary_to_superop
-
+from pygsti import SpaceT
 
 class GaugeGroup(_NicelySerializable):
     """
@@ -456,7 +456,7 @@ def transform_matrix(self):
         -------
         numpy.ndarray
         """
-        return self._operation.to_dense(on_space='minimal')
+        return self._operation.to_dense("minimal")
 
     @property
     def transform_matrix_inverse(self):
@@ -468,7 +468,7 @@ def transform_matrix_inverse(self):
         numpy.ndarray
         """
         if self._inv_matrix is None:
-            self._inv_matrix = _np.linalg.inv(self._operation.to_dense(on_space='minimal'))
+            self._inv_matrix = _np.linalg.inv(self._operation.to_dense("minimal"))
         return self._inv_matrix
 
     def deriv_wrt_params(self, wrt_filter=None):
diff --git a/pygsti/models/modelconstruction.py b/pygsti/models/modelconstruction.py
index df4cfc879..42c4763d2 100644
--- a/pygsti/models/modelconstruction.py
+++ b/pygsti/models/modelconstruction.py
@@ -47,6 +47,7 @@
 from pygsti.tools import internalgates as _itgs
 from pygsti.tools import optools as _ot
 from pygsti.tools import listtools as _lt
+from pygsti import SpaceT
 from pygsti.baseobjs.basisconstructors import sqrt2, id2x2, sigmax, sigmay, sigmaz
 from pygsti.baseobjs.verbosityprinter import VerbosityPrinter as _VerbosityPrinter
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
@@ -252,7 +253,7 @@ def to_labels(lbls):
             # a complex 2*num_qubits x 2*num_qubits mx unitary on full space in Pauli-product basis
             Uop_embed = _op.EmbeddedOp(state_space, labels, Uop)
             # a real 4*num_qubits x 4*num_qubits mx superoperator in final basis
-            superop_mx_pp = Uop_embed.to_dense(on_space='HilbertSchmidt')
+            superop_mx_pp = Uop_embed.to_dense("HilbertSchmidt")
             # a real 4*num_qubits x 4*num_qubits mx superoperator in final basis
             superop_mx_in_basis = _bt.change_basis(superop_mx_pp, 'pp', basis)
 
@@ -299,7 +300,7 @@ def to_labels(lbls):
             # a complex 2*num_qubits x 2*num_qubits mx unitary on full space in Pauli-product basis
             Uop_embed = _op.EmbeddedOp(state_space, (label,), Uop)
             # a real 4*num_qubits x 4*num_qubits mx superoperator in Pauli-product basis
-            superop_mx_pp = Uop_embed.to_dense(on_space='HilbertSchmidt')
+            superop_mx_pp = Uop_embed.to_dense("HilbertSchmidt")
             # a real 4*num_qubits x 4*num_qubits mx superoperator in final basis
             superop_mx_in_basis = _bt.change_basis(superop_mx_pp, 'pp', basis)
 
@@ -319,7 +320,7 @@ def to_labels(lbls):
             # a complex 2*num_qubits x 2*num_qubits mx unitary on full space in Pauli-product basis
             Uop_embed = _op.EmbeddedOp(state_space, (label,), Uop)
             # a real 4*num_qubits x 4*num_qubits mx superoperator in Pauli-product basis
-            superop_mx_pp = Uop_embed.to_dense(on_space='HilbertSchmidt')
+            superop_mx_pp = Uop_embed.to_dense("HilbertSchmidt")
             # a real 4*num_qubits x 4*num_qubits mx superoperator in final basis
             superop_mx_in_basis = _bt.change_basis(superop_mx_pp, 'pp', basis)
 
@@ -362,7 +363,7 @@ def to_labels(lbls):
             # a complex 2*num_qubits x 2*num_qubits mx unitary on full space
             Uop_embed = _op.EmbeddedOp(state_space, [label1, label2], Uop)
             # a real 4*num_qubits x 4*num_qubits mx superoperator in Pauli-product basis
-            superop_mx_pp = Uop_embed.to_dense(on_space='HilbertSchmidt')
+            superop_mx_pp = Uop_embed.to_dense("HilbertSchmidt")
             # a real 4*num_qubits x 4*num_qubits mx superoperator in final basis
             superop_mx_in_basis = _bt.change_basis(superop_mx_pp, 'pp', basis)
 
@@ -891,7 +892,7 @@ def _embed_unitary(statespace, target_labels, unitary):
                         raise NotImplementedError("'Iz' instrument can only be constructed on a space of *qubits*")
                     for ekey, effect_vec in _povm.ComputationalBasisPOVM(nqubits=len(qudit_labels), evotype=evotype,
                                                                          state_space=state_space).items():
-                        E = effect_vec.to_dense('HilbertSchmidt').reshape((state_space.dim, 1))
+                        E = effect_vec.to_dense("HilbertSchmidt").reshape((state_space.dim, 1))
                         inst_members[ekey] = _np.dot(E, E.T)  # (effect vector is a column vector)
                     ideal_instrument = _instrument.Instrument(inst_members)
                 else:
diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py
index 62266315f..a88975f16 100644
--- a/pygsti/objectivefns/objectivefns.py
+++ b/pygsti/objectivefns/objectivefns.py
@@ -27,7 +27,7 @@
 from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 from pygsti.baseobjs.verbosityprinter import VerbosityPrinter as _VerbosityPrinter
 from pygsti.models.model import OpModel as _OpModel
-
+from pygsti import SpaceT
 
 def _objfn(objfn_cls, model, dataset, circuits=None,
            regularization=None, penalties=None, op_label_aliases=None,
@@ -5790,11 +5790,11 @@ def _spam_penalty(mdl, prefactor, op_basis):
     return prefactor * (_np.sqrt(
         _np.array([
             _tools.tracenorm(
-                _tools.vec_to_stdmx(prepvec.to_dense(on_space='minimal'), op_basis)
+                _tools.vec_to_stdmx(prepvec.to_dense("minimal"), op_basis)
             ) for prepvec in mdl.preps.values()
         ] + [
             _tools.tracenorm(
-                _tools.vec_to_stdmx(mdl.povms[plbl][elbl].to_dense(on_space='minimal'), op_basis)
+                _tools.vec_to_stdmx(mdl.povms[plbl][elbl].to_dense("minimal"), op_basis)
             ) for plbl in mdl.povms for elbl in mdl.povms[plbl]], 'd')
     ))
 
diff --git a/pygsti/pgtypes.py b/pygsti/pgtypes.py
new file mode 100644
index 000000000..5de0cf0fe
--- /dev/null
+++ b/pygsti/pgtypes.py
@@ -0,0 +1,3 @@
+from typing import Literal
+
+SpaceT = Literal["minimal", "Hilbert", "HilbertSchmidt"]
\ No newline at end of file
diff --git a/pygsti/protocols/freeformsim.py b/pygsti/protocols/freeformsim.py
index 90b15bc3a..552cbd527 100644
--- a/pygsti/protocols/freeformsim.py
+++ b/pygsti/protocols/freeformsim.py
@@ -16,7 +16,7 @@
 from pygsti.circuits.circuit import Circuit as _Circuit
 from pygsti.data.freedataset import FreeformDataSet as _FreeformDataSet
 from pygsti.modelmembers import states as _state
-
+from pygsti import SpaceT
 
 class FreeformDataSimulator(_proto.DataSimulator):
     """
@@ -114,7 +114,7 @@ def compute_process_matrix(self, model, circuit, include_final_state=False, incl
         if include_final_state or include_probabilities:
             ret = [mx]
             rho = model.circuit_layer_operator(prep, 'prep')
-            final_state = _state.StaticState(_np.dot(mx, rho.to_dense(on_space='HilbertSchmidt')),
+            final_state = _state.StaticState(_np.dot(mx, rho.to_dense("HilbertSchmidt")),
                                              model.basis, model.evotype, model.state_space)
             if include_final_state:
                 ret.append(final_state)
diff --git a/pygsti/report/modelfunction.py b/pygsti/report/modelfunction.py
index a4c0d8c2f..6a2d89da2 100644
--- a/pygsti/report/modelfunction.py
+++ b/pygsti/report/modelfunction.py
@@ -14,7 +14,7 @@
 from pygsti.models.localnoisemodel import LocalNoiseModel as _LocalNoiseModel
 from pygsti.baseobjs.basis import Basis as _Basis
 from pygsti.baseobjs.basis import TensorProdBasis as _TensorProdBasis
-
+from pygsti import SpaceT
 
 class ModelFunction(object):
     """
@@ -195,7 +195,7 @@ def __init__(self, model, gl, *args, **kwargs):
 
         def evaluate(self, model):
             """ Evaluate this gate-set-function at `model`."""
-            return fn(model.operations[self.gl].to_dense(on_space='HilbertSchmidt'), model.basis,
+            return fn(model.operations[self.gl].to_dense("HilbertSchmidt"), model.basis,
                       *self.args, **self.kwargs)
 
     GSFTemp.__name__ = fn.__name__ + str("_class")
@@ -239,8 +239,8 @@ def __init__(self, model1, model2, gl, *args, **kwargs):
         def evaluate(self, model):
             """ Evaluate this gate-set-function at `model`."""
             if isinstance(model, _ExplicitOpModel):
-                return fn(model.operations[self.gl].to_dense(on_space='HilbertSchmidt'),
-                          self.other_model.operations[self.gl].to_dense(on_space='HilbertSchmidt'),
+                return fn(model.operations[self.gl].to_dense("HilbertSchmidt"),
+                          self.other_model.operations[self.gl].to_dense("HilbertSchmidt"),
                           model.basis, *self.args, **self.kwargs)  # assume functions want *dense* gates
             elif isinstance(model, _LocalNoiseModel):
                 opdim = model.operation_blks['gates'][self.gl].dim
@@ -253,8 +253,8 @@ def evaluate(self, model):
                 else:
                     raise ValueError(f"Could not convert model basis (name={model.basis.name}) to an appropriate {opdim}-dim basis!")
 
-                return fn(model.operation_blks['gates'][self.gl].to_dense(on_space='HilbertSchmidt'),
-                          self.other_model.operation_blks['gates'][self.gl].to_dense(on_space='HilbertSchmidt'),
+                return fn(model.operation_blks['gates'][self.gl].to_dense("HilbertSchmidt"),
+                          self.other_model.operation_blks['gates'][self.gl].to_dense("HilbertSchmidt"),
                           basis, *self.args, **self.kwargs)  # assume functions want *dense* gates
             else:
                 raise ValueError(f"Unsupported model type: {type(model)}!")
@@ -344,11 +344,11 @@ def __init__(self, model, lbl, typ, *args, **kwargs):
         def evaluate(self, model):
             """ Evaluate this gate-set-function at `model`."""
             if self.typ == "prep":
-                return fn(model.preps[self.lbl].to_dense(on_space='HilbertSchmidt'), model.basis,
+                return fn(model.preps[self.lbl].to_dense("HilbertSchmidt"), model.basis,
                           *self.args, **self.kwargs)
             else:
                 povmlbl, Elbl = self.lbl.split(":")  # for effect, lbl must == "povmLbl:ELbl"
-                return fn(model.povms[povmlbl][Elbl].to_dense(on_space='HilbertSchmidt'), model.basis,
+                return fn(model.povms[povmlbl][Elbl].to_dense("HilbertSchmidt"), model.basis,
                           *self.args, **self.kwargs)
 
     GSFTemp.__name__ = fn.__name__ + str("_class")
@@ -399,13 +399,13 @@ def __init__(self, model1, model2, lbl, typ, *args, **kwargs):
         def evaluate(self, model):
             """ Evaluate this gate-set-function at `model`."""
             if self.typ == "prep":
-                return fn(model.preps[self.lbl].to_dense(on_space='HilbertSchmidt'),
-                          self.other_vecsrc[self.lbl].to_dense(on_space='HilbertSchmidt'),
+                return fn(model.preps[self.lbl].to_dense("HilbertSchmidt"),
+                          self.other_vecsrc[self.lbl].to_dense("HilbertSchmidt"),
                           model.basis, *self.args, **self.kwargs)
             else:
                 povmlbl, Elbl = self.lbl.split(":")  # for effect, lbl must == "povmLbl:ELbl"
-                return fn(model.povms[povmlbl][Elbl].to_dense(on_space='HilbertSchmidt'),
-                          self.other_vecsrc[povmlbl][Elbl].to_dense(on_space='HilbertSchmidt'),
+                return fn(model.povms[povmlbl][Elbl].to_dense("HilbertSchmidt"),
+                          self.other_vecsrc[povmlbl][Elbl].to_dense("HilbertSchmidt"),
                           model.basis, *self.args, **self.kwargs)
 
     GSFTemp.__name__ = fn.__name__ + str("_class")
diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py
index dfa46a181..cd1b403ee 100644
--- a/pygsti/report/reportables.py
+++ b/pygsti/report/reportables.py
@@ -30,7 +30,7 @@
 from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as _LEEL
 from pygsti.modelmembers.operations.lindbladcoefficients import LindbladCoefficientBlock as _LindbladCoefficientBlock
 from pygsti.models.explicitmodel import ExplicitOpModel as _ExplicitOpModel
-
+from pygsti import SpaceT
 
 _CVXPY_AVAILABLE = importlib.util.find_spec('cvxpy') is not None
 
@@ -121,8 +121,8 @@ def spam_dotprods(rho_vecs, povms):
         j = 0
         for povm in povms:
             for EVec in povm.values():
-                ret[i, j] = _np.vdot(EVec.to_dense(on_space='HilbertSchmidt'),
-                                     rhoVec.to_dense(on_space='HilbertSchmidt'))
+                ret[i, j] = _np.vdot(EVec.to_dense("HilbertSchmidt"),
+                                     rhoVec.to_dense("HilbertSchmidt"))
                 j += 1
                 # to_dense() gives a 1D array, so no need to transpose EVec
     return ret
@@ -230,7 +230,7 @@ def evaluate(self, model):
         -------
         numpy.ndarray
         """
-        evals, evecs = _np.linalg.eig(model.operations[self.oplabel].to_dense(on_space='HilbertSchmidt'))
+        evals, evecs = _np.linalg.eig(model.operations[self.oplabel].to_dense("HilbertSchmidt"))
 
         ev_list = list(enumerate(evals))
         ev_list.sort(key=lambda tup: abs(tup[1]), reverse=True)
@@ -1031,7 +1031,7 @@ def angles_btwn_rotn_axes(model):
     angles_btwn_rotn_axes = _np.zeros((len(opLabels), len(opLabels)), 'd')
 
     for i, gl in enumerate(opLabels):
-        decomp = _tools.decompose_gate_matrix(model.operations[gl].to_dense(on_space='HilbertSchmidt'))
+        decomp = _tools.decompose_gate_matrix(model.operations[gl].to_dense("HilbertSchmidt"))
         rotnAngle = decomp.get('pi rotations', 'X')
         axisOfRotn = decomp.get('axis of rotation', None)
 
@@ -1222,9 +1222,9 @@ class HalfDiamondNorm(_modf.ModelFunction):
         def __init__(self, model_a, model_b, oplabel):
             self.oplabel = oplabel
             if isinstance(model_b, _ExplicitOpModel):
-                self.B = model_b.operations[oplabel].to_dense(on_space='HilbertSchmidt')
+                self.B = model_b.operations[oplabel].to_dense("HilbertSchmidt")
             else:
-                self.B = model_b.operation_blks['gates'][oplabel].to_dense(on_space='HilbertSchmidt')
+                self.B = model_b.operation_blks['gates'][oplabel].to_dense("HilbertSchmidt")
             self.d = int(round(_np.sqrt(model_a.dim)))
             _modf.ModelFunction.__init__(self, model_a, [("gate", oplabel)])
 
@@ -1243,10 +1243,10 @@ def evaluate(self, model):
             """
             gl = self.oplabel
             if isinstance(model, _ExplicitOpModel):
-                dm, W = _tools.diamonddist(model.operations[gl].to_dense(on_space='HilbertSchmidt'),
+                dm, W = _tools.diamonddist(model.operations[gl].to_dense("HilbertSchmidt"),
                                            self.B, model.basis, return_x=True)
             else:
-                dm, W = _tools.diamonddist(model.operation_blks['gates'][gl].to_dense(on_space='HilbertSchmidt'),
+                dm, W = _tools.diamonddist(model.operation_blks['gates'][gl].to_dense("HilbertSchmidt"),
                                            self.B, 'pp', return_x=True)  # HACK - need to get basis from model 'pp' HARDCODED for now
 
             self.W = W
@@ -1267,9 +1267,9 @@ def evaluate_nearby(self, nearby_model):
             """
             mxBasis = nearby_model.basis
             if isinstance(nearby_model, _ExplicitOpModel):
-                A = nearby_model.operations[self.oplabel].to_dense(on_space='HilbertSchmidt')
+                A = nearby_model.operations[self.oplabel].to_dense("HilbertSchmidt")
             else:
-                A = nearby_model.operation_blks['gates'][self.oplabel].to_dense(on_space='HilbertSchmidt')
+                A = nearby_model.operation_blks['gates'][self.oplabel].to_dense("HilbertSchmidt")
                 mxBasis = 'pp'  # HACK need to set mxBasis based on model but not the full model basis
             JAstd = self.d * _tools.fast_jamiolkowski_iso_std(A, mxBasis)
             JBstd = self.d * _tools.fast_jamiolkowski_iso_std(self.B, mxBasis)
@@ -2146,8 +2146,8 @@ def general_decomposition(model_a, model_b):
     mxBasis = model_b.basis  # B is usually the target which has a well-defined basis
 
     for gl in opLabels:
-        gate = model_a.operations[gl].to_dense(on_space='HilbertSchmidt')
-        targetOp = model_b.operations[gl].to_dense(on_space='HilbertSchmidt')
+        gate = model_a.operations[gl].to_dense("HilbertSchmidt")
+        targetOp = model_b.operations[gl].to_dense("HilbertSchmidt")
         gl = str(gl)  # Label -> str for decomp-dict keys
 
         target_evals = _np.linalg.eigvals(targetOp)
@@ -2621,8 +2621,8 @@ def instrument_half_diamond_norm(a, b, mx_basis):
         aa, bb = i * adim, (i + 1) * adim
         for j in range(nComps):
             cc, dd = j * adim, (j + 1) * adim
-            composite_op[aa:bb, cc:dd] = a[clbl].to_dense(on_space='HilbertSchmidt')
-            composite_top[aa:bb, cc:dd] = b[clbl].to_dense(on_space='HilbertSchmidt')
+            composite_op[aa:bb, cc:dd] = a[clbl].to_dense("HilbertSchmidt")
+            composite_top[aa:bb, cc:dd] = b[clbl].to_dense("HilbertSchmidt")
     return half_diamond_norm(composite_op, composite_top, sumbasis)
 
 
diff --git a/pygsti/report/workspacetables.py b/pygsti/report/workspacetables.py
index 9c358b592..4147af167 100644
--- a/pygsti/report/workspacetables.py
+++ b/pygsti/report/workspacetables.py
@@ -33,6 +33,7 @@
 from pygsti.circuits.circuit import Circuit as _Circuit
 from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as _LEEL
 from pygsti.data import DataSet as _DataSet
+from pygsti import SpaceT
 
 
 class BlankTable(WorkspaceTable):
@@ -528,7 +529,7 @@ def _create(self, models, titles, display_as, confidence_region_infos):
                     intervalMx = None
 
                 #turn the op matrix into a ReportableQty
-                op_matrix = _ReportableQty(op.to_dense('HilbertSchmidt'), errbar=intervalMx)
+                op_matrix = _ReportableQty(op.to_dense("HilbertSchmidt"), errbar=intervalMx)
 
                 if display_as == "numbers":
                     row_data.append(op_matrix)
@@ -835,8 +836,8 @@ def _get_gig_decomp(mx, tmx):  # "Gauge invariant gateset" decomposition
         op_decomps = {}
         for gl in opLabels:
             try:
-                op_decomps[gl] = _get_gig_decomp(model.operations[gl].to_dense(on_space='HilbertSchmidt'),
-                                                 target_model.operations[gl].to_dense(on_space='HilbertSchmidt'))
+                op_decomps[gl] = _get_gig_decomp(model.operations[gl].to_dense("HilbertSchmidt"),
+                                                 target_model.operations[gl].to_dense("HilbertSchmidt"))
                 M = max(M, max(_np.abs((op_decomps[gl][1] - I).flat)))  # update max
             except Exception as e:
                 _warnings.warn("Failed gauge-robust decomposition of %s op:\n%s" % (gl, str(e)))
@@ -1018,8 +1019,8 @@ def _create(self, model, target_model, metric, confidence_region_info):
         mdl_in_best_gauge = []
         target_mdl_in_best_gauge = []
         for lbl in opLabels:
-            gate_mx = orig_model.operations[lbl].to_dense(on_space='HilbertSchmidt')
-            target_gate_mx = target_model.operations[lbl].to_dense(on_space='HilbertSchmidt')
+            gate_mx = orig_model.operations[lbl].to_dense("HilbertSchmidt")
+            target_gate_mx = target_model.operations[lbl].to_dense("HilbertSchmidt")
             Ugauge = _tools.compute_best_case_gauge_transform(gate_mx, target_gate_mx)
             Ugg = _models.gaugegroup.FullGaugeGroupElement(_np.linalg.inv(Ugauge))
             # transforms gates as Ugauge * gate * Ugauge_inv
@@ -2557,7 +2558,7 @@ def _create(self, model, target_model,
 
                 if isinstance(gl, _baseobjs.Label) or isinstance(gl, str):
                     # no error bars
-                    target_evals = _np.linalg.eigvals(target_model.operations[gl].to_dense(on_space='HilbertSchmidt'))
+                    target_evals = _np.linalg.eigvals(target_model.operations[gl].to_dense("HilbertSchmidt"))
                 else:
                     target_evals = _np.linalg.eigvals(target_model.sim.product(gl))  # no error bars
 
@@ -2672,11 +2673,11 @@ def _create(self, model, target_model,
                 row_formatters = [None]
 
                 #FUTURE: use reportables to get instrument eigenvalues
-                evals = _ReportableQty(_np.linalg.eigvals(comp.to_dense(on_space='HilbertSchmidt')))
+                evals = _ReportableQty(_np.linalg.eigvals(comp.to_dense("HilbertSchmidt")))
                 evals = evals.reshape(evals.size, 1)
 
                 if target_model is not None:
-                    target_evals = _np.linalg.eigvals(tcomp.to_dense(on_space='HilbertSchmidt'))  # no error bars
+                    target_evals = _np.linalg.eigvals(tcomp.to_dense("HilbertSchmidt"))  # no error bars
                     #Note: no support for relative eigenvalues of instruments (yet)
 
                     # permute target eigenvalues according to min-weight matching
diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py
index 3cea8291d..c53ec9d7a 100644
--- a/pygsti/tools/jamiolkowski.py
+++ b/pygsti/tools/jamiolkowski.py
@@ -15,7 +15,7 @@
 from pygsti.tools import basistools as _bt
 from pygsti.tools import matrixtools as _mt
 from pygsti.baseobjs.basis import Basis as _Basis
-
+from pygsti import SpaceT
 
 # Gate Mx G:      rho  --> G rho                    where G and rho are in the Pauli basis (by definition/convention)                                           # noqa
 #            vec(rhoS) --> GStd vec(rhoS)           where GS and rhoS are in the std basis, GS = PtoS * G * StoP                                                # noqa
@@ -348,7 +348,7 @@ def sums_of_negative_choi_eigenvalues(model):
     """
     ret = []
     for (_, gate) in model.operations.items():
-        J = fast_jamiolkowski_iso_std(gate.to_dense(on_space='HilbertSchmidt'), model.basis)  # Choi mx basis doesn't matter
+        J = fast_jamiolkowski_iso_std(gate.to_dense("HilbertSchmidt"), model.basis)  # Choi mx basis doesn't matter
         evals = _np.linalg.eigvals(J)  # could use eigvalsh, but wary of this since eigh can be wrong...
         sumOfNeg = 0.0
         for ev in evals:
diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py
index 598a72605..a4f2feec1 100644
--- a/pygsti/tools/optools.py
+++ b/pygsti/tools/optools.py
@@ -27,6 +27,7 @@
 from pygsti.baseobjs.label import Label as _Label
 from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
+from pygsti import SpaceT
 
 IMAG_TOL = 1e-7  # tolerance for imaginary part being considered zero
 
@@ -1087,8 +1088,8 @@ def instrument_diamonddist(a, b, mx_basis):
         aa, bb = i * adim, (i + 1) * adim
         for j in range(nComps):
             cc, dd = j * adim, (j + 1) * adim
-            composite_op[aa:bb, cc:dd] = a[clbl].to_dense(on_space='HilbertSchmidt')
-            composite_top[aa:bb, cc:dd] = b[clbl].to_dense(on_space='HilbertSchmidt')
+            composite_op[aa:bb, cc:dd] = a[clbl].to_dense("HilbertSchmidt")
+            composite_top[aa:bb, cc:dd] = b[clbl].to_dense("HilbertSchmidt")
     return diamonddist(composite_op, composite_top, sumbasis)
 
 
@@ -2654,8 +2655,8 @@ def project_to_target_eigenspace(model, target_model, eps=1e-6):
         #Essentially, we want to replace the eigenvalues of `tgt_gate`
         # (and *only* the eigenvalues) with those of `gate`.  This is what
         # a "best gate gauge transform does" (by definition)
-        gate_mx = gate.to_dense(on_space='minimal')
-        Ugauge = compute_best_case_gauge_transform(gate_mx, tgt_gate.to_dense(on_space='minimal'))
+        gate_mx = gate.to_dense("minimal")
+        Ugauge = compute_best_case_gauge_transform(gate_mx, tgt_gate.to_dense("minimal"))
         Ugauge_inv = _np.linalg.inv(Ugauge)
 
         epgate = _np.dot(Ugauge, _np.dot(gate_mx, Ugauge_inv))
diff --git a/pygsti/tools/rbtheory.py b/pygsti/tools/rbtheory.py
index 9649fc6c4..bb0cd03ed 100644
--- a/pygsti/tools/rbtheory.py
+++ b/pygsti/tools/rbtheory.py
@@ -17,7 +17,7 @@
 from pygsti.tools import matrixtools as _mtls
 from pygsti.tools import optools as _optls
 from pygsti.tools import rbtools as _rbtls
-
+from pygsti import SpaceT
 
 def predicted_rb_number(model, target_model, weights=None, d=None, rtype='EI'):
     """
@@ -328,8 +328,8 @@ def L_matrix(model, target_model, weights=None):  # noqa N802
     normalizer = _np.sum(_np.array([weights[key] for key in list(target_model.operations.keys())]))
     L_matrix = (1 / normalizer) * _np.sum(
         weights[key] * _np.kron(
-            model.operations[key].to_dense(on_space='HilbertSchmidt').T,
-            _np.linalg.inv(target_model.operations[key].to_dense(on_space='HilbertSchmidt'))
+            model.operations[key].to_dense("HilbertSchmidt").T,
+            _np.linalg.inv(target_model.operations[key].to_dense("HilbertSchmidt"))
         ) for key in target_model.operations.keys())
 
     return L_matrix
diff --git a/test/unit/modelmembers/test_kraus_interface.py b/test/unit/modelmembers/test_kraus_interface.py
index 625e55df7..ab0b27fd1 100644
--- a/test/unit/modelmembers/test_kraus_interface.py
+++ b/test/unit/modelmembers/test_kraus_interface.py
@@ -171,13 +171,13 @@ def setUp(self):
                                                [0., 0.9, 0., 0.],
                                                [0., 0., 0.9, 0.],
                                                [0., 0., 0., 0.9]])
-        self.assertArraysAlmostEqual(op.to_dense(on_space='HilbertSchmidt'), self.expected_idle_superop)
+        self.assertArraysAlmostEqual(op.to_dense('HilbertSchmidt'), self.expected_idle_superop)
 
     def test_stochastic_op_creation(self):
         ss = QubitSpace(1)
         op = StochasticNoiseOp(ss, initial_rates=[0.025, 0.025, 0.025], evotype=self.evotype) # 0.025 = 0.1/4
         try:
-            self.assertArraysAlmostEqual(op.to_dense(on_space='HilbertSchmidt'), self.expected_idle_superop)
+            self.assertArraysAlmostEqual(op.to_dense('HilbertSchmidt'), self.expected_idle_superop)
         except NotImplementedError:
             pass  # ok if to_dense not implemented, as for CHP evotype
 
@@ -198,7 +198,7 @@ def test_depol_model(self):
         self.assertTrue(isinstance(ops[('Gypi2', 0)], ComposedOp))
 
         try:
-            Gx_error = mdl_sto.operations['Gxpi2', 0].factorops[1].to_dense(on_space='HilbertSchmidt')
+            Gx_error = mdl_sto.operations['Gxpi2', 0].factorops[1].to_dense('HilbertSchmidt')
             self.assertArraysAlmostEqual(Gx_error, self.expected_idle_superop)
         except NotImplementedError:
             pass  # ok if not implemented, as for CHP evotype