"""
Defines classes with represent SPAM operations, along with supporting
functionality.
"""
#***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License.  You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************

import numpy as _np
import scipy.sparse as _sps
import collections as _collections
import numbers as _numbers
import itertools as _itertools
import functools as _functools
import copy as _copy

from .. import optimize as _opt
from ..tools import matrixtools as _mt
from ..tools import optools as _gt
from ..tools import basistools as _bt
from ..tools import listtools as _lt
from ..tools import slicetools as _slct
from ..tools import compattools as _compat
from ..tools import symplectic as _symp
from .basis import Basis as _Basis
from .protectedarray import ProtectedArray as _ProtectedArray
from . import modelmember as _modelmember

from . import term as _term
from . import stabilizer as _stabilizer
from .polynomial import Polynomial as _Polynomial
from . import replib
from .opcalc import bulk_eval_compact_polys_complex as _bulk_eval_compact_polys_complex

IMAG_TOL = 1e-8  # tolerance for imaginary part being considered zero


def optimize_spamvec(vecToOptimize, targetVec):
    """
    Optimize the parameters of vecToOptimize so that the
      the resulting SPAM vector is as close as possible to
      targetVec.

    This is trivial for the case of FullSPAMVec
      instances, but for other types of parameterization
      this involves an iterative optimization over all the
      parameters of vecToOptimize.

    Parameters
    ----------
    vecToOptimize : SPAMVec
      The vector to optimize. This object gets altered.

    targetVec : SPAMVec
      The SPAM vector used as the target.

    Returns
    -------
    None
    """

    #TODO: cleanup this code:
    if isinstance(vecToOptimize, StaticSPAMVec):
        return  # nothing to optimize

    if isinstance(vecToOptimize, FullSPAMVec):
        if(targetVec.dim != vecToOptimize.dim):  # special case: gates can have different overall dimension
            vecToOptimize.dim = targetVec.dim  # this is a HACK to allow model selection code to work correctly
        vecToOptimize.set_value(targetVec)  # just copy entire overall matrix since fully parameterized
        return

    assert(targetVec.dim == vecToOptimize.dim)  # vectors must have the same overall dimension
    targetVector = _np.asarray(targetVec)

    def _objective_func(param_vec):
        vecToOptimize.from_vector(param_vec)
        return _mt.frobeniusnorm(vecToOptimize - targetVector)

    x0 = vecToOptimize.to_vector()
    minSol = _opt.minimize(_objective_func, x0, method='BFGS', maxiter=10000, maxfev=10000,
                           tol=1e-6, callback=None)

    vecToOptimize.from_vector(minSol.x)
    #print("DEBUG: optimized vector to min frobenius distance %g" % _mt.frobeniusnorm(vecToOptimize-targetVector))


def convert(spamvec, toType, basis, extra=None):
    """
    Convert SPAM vector to a new type of parameterization, potentially
    creating a new SPAMVec object.  Raises ValueError for invalid conversions.

    Parameters
    ----------
    spamvec : SPAMVec
        SPAM vector to convert

    toType : {"full","TP","static","static unitary","clifford",LINDBLAD}
        The type of parameterizaton to convert to.  "LINDBLAD" is a placeholder
        for the various Lindblad parameterization types.  See
        :method:`Model.set_all_parameterizations` for more details.

    basis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The basis for `spamvec`.  Allowed values are Matrix-unit (std),
        Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt)
        (or a custom basis object).

    extra : object, optional
        Additional information for conversion.


    Returns
    -------
    SPAMVec
       The converted SPAM vector, usually a distinct
       object from the object passed as input.
    """
    if toType == "full":
        if isinstance(spamvec, FullSPAMVec):
            return spamvec  # no conversion necessary
        else:
            typ = spamvec._prep_or_effect if isinstance(spamvec, SPAMVec) else "prep"
            return FullSPAMVec(spamvec.todense(), typ=typ)

    elif toType == "TP":
        if isinstance(spamvec, TPSPAMVec):
            return spamvec  # no conversion necessary
        else:
            return TPSPAMVec(spamvec.todense())
            # above will raise ValueError if conversion cannot be done

    elif toType == "TrueCPTP":  # a non-lindbladian CPTP spamvec that hasn't worked well...
        if isinstance(spamvec, CPTPSPAMVec):
            return spamvec  # no conversion necessary
        else:
            return CPTPSPAMVec(spamvec, basis)
            # above will raise ValueError if conversion cannot be done

    elif toType == "static":
        if isinstance(spamvec, StaticSPAMVec):
            return spamvec  # no conversion necessary
        else:
            typ = spamvec._prep_or_effect if isinstance(spamvec, SPAMVec) else "prep"
            return StaticSPAMVec(spamvec, typ=typ)

    elif toType == "static unitary":
        dmvec = _bt.change_basis(spamvec.todense(), basis, 'std')
        purevec = _gt.dmvec_to_state(dmvec)
        return StaticSPAMVec(purevec, "statevec", spamvec._prep_or_effect)

    elif _gt.is_valid_lindblad_paramtype(toType):

        if extra is None:
            purevec = spamvec  # right now, we don't try to extract a "closest pure vec"
            # to spamvec - below will fail if spamvec isn't pure.
        else:
            purevec = extra  # assume extra info is a pure vector

        nQubits = _np.log2(spamvec.dim) / 2.0
        bQubits = bool(abs(nQubits - round(nQubits)) < 1e-10)  # integer # of qubits?
        proj_basis = "pp" if (basis == "pp" or bQubits) else basis
        typ = spamvec._prep_or_effect if isinstance(spamvec, SPAMVec) else "prep"

        return LindbladSPAMVec.from_spamvec_obj(
            spamvec, typ, toType, None, proj_basis, basis,
            truncate=True, lazy=True)

    elif toType == "clifford":
        if isinstance(spamvec, StabilizerSPAMVec):
            return spamvec  # no conversion necessary

        purevec = spamvec.flatten()  # assume a pure state (otherwise would
        # need to change Model dim)
        return StabilizerSPAMVec.from_dense_purevec(purevec)

    else:
        raise ValueError("Invalid toType argument: %s" % toType)


def _convert_to_lindblad_base(vec, typ, new_evotype, mxBasis="pp"):
    """
    Attempts to convert `vec` to a static (0 params) SPAMVec with
    evoution type `new_evotype`.  Used to convert spam vecs to
    being LindbladSPAMVec objects.
    """
    if vec._evotype == new_evotype and vec.num_params() == 0:
        return vec  # no conversion necessary
    if new_evotype == "densitymx":
        return StaticSPAMVec(vec.todense(), "densitymx", typ)
    if new_evotype in ("svterm", "cterm"):
        if isinstance(vec, ComputationalSPAMVec):  # special case when conversion is easy
            return ComputationalSPAMVec(vec._zvals, new_evotype, typ)
        elif vec._evotype == "densitymx":
            # then try to extract a (static) pure state from vec wth
            # evotype 'statevec' or 'stabilizer' <=> 'svterm', 'cterm'
            if isinstance(vec, DenseSPAMVec):
                dmvec = _bt.change_basis(vec, mxBasis, 'std')
                purestate = StaticSPAMVec(_gt.dmvec_to_state(dmvec), 'statevec', typ)
            elif isinstance(vec, PureStateSPAMVec):
                purestate = vec.pure_state_vec  # evotype 'statevec'
            else:
                raise ValueError("Unable to obtain pure state from density matrix type %s!" % type(vec))

            if new_evotype == "cterm":  # then purestate 'statevec' => 'stabilizer' (if possible)
                if typ == "prep":
                    purestate = StabilizerSPAMVec.from_dense_purevec(purestate.todense())
                else:  # type == "effect"
                    purestate = StabilizerEffectVec.from_dense_purevec(purestate.todense())

            return PureStateSPAMVec(purestate, new_evotype, mxBasis, typ)

    raise ValueError("Could not convert %s (evotype %s) to %s w/0 params!" %
                     (str(type(vec)), vec._evotype, new_evotype))


def finite_difference_deriv_wrt_params(spamvec, wrtFilter=None, eps=1e-7):
    """
    Computes a finite-difference Jacobian for a SPAMVec object.

    The returned value is a matrix whose columns are the vectorized
    derivatives of the spam vector with respect to a single
    parameter, matching the format expected from the spam vectors's
    `deriv_wrt_params` method.

    Parameters
    ----------
    spamvec : SPAMVec
        The spam vector object to compute a Jacobian for.

    eps : float, optional
        The finite difference step to use.

    Returns
    -------
    numpy.ndarray
        An M by N matrix where M is the number of gate elements and
        N is the number of gate parameters.
    """
    dim = spamvec.get_dimension()
    spamvec2 = spamvec.copy()
    p = spamvec.to_vector()
    fd_deriv = _np.empty((dim, spamvec.num_params()), 'd')  # assume real (?)

    for i in range(spamvec.num_params()):
        p_plus_dp = p.copy()
        p_plus_dp[i] += eps
        spamvec2.from_vector(p_plus_dp, close=True)
        fd_deriv[:, i:i + 1] = (spamvec2 - spamvec) / eps

    fd_deriv.shape = [dim, spamvec.num_params()]
    if wrtFilter is None:
        return fd_deriv
    else:
        return _np.take(fd_deriv, wrtFilter, axis=1)


def check_deriv_wrt_params(spamvec, deriv_to_check=None, wrtFilter=None, eps=1e-7):
    """
    Checks the `deriv_wrt_params` method of a SPAMVec object.

    This routine is meant to be used as an aid in testing and debugging
    SPAMVec classes by comparing the finite-difference Jacobian that
    *should* be returned by `spamvec.deriv_wrt_params` with the one that
    actually is.  A ValueError is raised if the two do not match.

    Parameters
    ----------
    spamvec : SPAMVec
        The gate object to test.

    deriv_to_check : numpy.ndarray or None, optional
        If not None, the Jacobian to compare against the finite difference
        result.  If None, `spamvec.deriv_wrt_parms()` is used.  Setting this
        argument can be useful when the function is called *within* a LinearOperator
        class's `deriv_wrt_params()` method itself as a part of testing.

    eps : float, optional
        The finite difference step to use.

    Returns
    -------
    None
    """
    fd_deriv = finite_difference_deriv_wrt_params(spamvec, wrtFilter, eps)
    if deriv_to_check is None:
        deriv_to_check = spamvec.deriv_wrt_params()

    #print("Deriv shapes = %s and %s" % (str(fd_deriv.shape),
    #                                    str(deriv_to_check.shape)))
    #print("finite difference deriv = \n",fd_deriv)
    #print("deriv_wrt_params deriv = \n",deriv_to_check)
    #print("deriv_wrt_params - finite diff deriv = \n",
    #      deriv_to_check - fd_deriv)

    for i in range(deriv_to_check.shape[0]):
        for j in range(deriv_to_check.shape[1]):
            diff = abs(deriv_to_check[i, j] - fd_deriv[i, j])
            if diff > 5 * eps:
                print("deriv_chk_mismatch: (%d,%d): %g (comp) - %g (fd) = %g" %
                      (i, j, deriv_to_check[i, j], fd_deriv[i, j], diff))

    if _np.linalg.norm(fd_deriv - deriv_to_check) > 100 * eps:
        raise ValueError("Failed check of deriv_wrt_params:\n"
                         " norm diff = %g" %
                         _np.linalg.norm(fd_deriv - deriv_to_check))


class SPAMVec(_modelmember.ModelMember):
    """
    Excapulates a parameterization of a state preparation OR POVM effect
    vector. This class is the  common base class for all specific
    parameterizations of a SPAM vector.
    """

    def __init__(self, rep, evotype, typ):
        """ Initialize a new SPAM Vector """
        if isinstance(rep, int):  # For operators that have no representation themselves (term ops)
            dim = rep             # allow passing an integer as `rep`.
            rep = None
        else:
            dim = rep.dim
        super(SPAMVec, self).__init__(dim, evotype)
        self._rep = rep
        self._prep_or_effect = typ

    @property
    def size(self):
        """
        Return the number of independent elements in this gate (when viewed as a dense array)
        """
        return self.dim

    @property
    def outcomes(self):
        """
        Return the z-value outcomes corresponding to this effect SPAM vector
        in the context of a stabilizer-state simulation.
        """
        raise NotImplementedError("'outcomes' property is not implemented for %s objects" % self.__class__.__name__)

    def set_value(self, vec):
        """
        Attempts to modify SPAMVec parameters so that the specified raw
        SPAM vector becomes vec.  Will raise ValueError if this operation
        is not possible.

        Parameters
        ----------
        vec : array_like or SPAMVec
            A numpy array representing a SPAM vector, or a SPAMVec object.

        Returns
        -------
        None
        """
        raise ValueError("Cannot set the value of a %s directly!" % self.__class__.__name__)

    def set_time(self, t):
        """
        Sets the current time for a time-dependent operator.  For time-independent
        operators (the default), this function does absolutely nothing.

        Parameters
        ----------
        t : float
            The current time.

        Returns
        -------
        None
        """
        pass

    def todense(self, scratch=None):
        """
        Return this SPAM vector as a (dense) numpy array.  The memory
        in `scratch` maybe used when it is not-None.
        """
        raise NotImplementedError("todense(...) not implemented for %s objects!" % self.__class__.__name__)

#    def torep(self, typ, outrep=None):
#        """
#        Return a "representation" object for this SPAM vector.
#        Such objects are primarily used internally by pyGSTi to compute
#        things like probabilities more efficiently.
#
#        Parameters
#        ----------
#        typ : {'prep','effect'}
#            The type of representation (for cases when the vector type is
#            not already defined).
#
#        outrep : StateRep
#            If not None, an existing state representation appropriate to this
#            SPAM vector that may be used instead of allocating a new one.
#
#        Returns
#        -------
#        StateRep
#        """
#        if typ == "prep":
#            if self._evotype == "statevec":
#                return replib.SVStateRep(self.todense())
#            elif self._evotype == "densitymx":
#                return replib.DMStateRep(self.todense())
#            raise NotImplementedError("torep(%s) not implemented for %s objects!" %
#                                      (self._evotype, self.__class__.__name__))
#        elif typ == "effect":
#            if self._evotype == "statevec":
#                return replib.SVEffectRep_Dense(self.todense())
#            elif self._evotype == "densitymx":
#                return replib.DMEffectRep_Dense(self.todense())
#            raise NotImplementedError("torep(%s) not implemented for %s objects!" %
#                                      (self._evotype, self.__class__.__name__))
#        else:
#            raise ValueError("Invalid `typ` argument for torep(): %s" % typ)

    def get_taylor_order_terms(self, order, max_poly_vars=100, return_poly_coeffs=False):
        """
        Get the `order`-th order Taylor-expansion terms of this SPAM vector.

        This function either constructs or returns a cached list of the terms at
        the given order.  Each term is "rank-1", meaning that it is a state
        preparation followed by or POVM effect preceded by actions on a
        density matrix `rho` of the form:

        `rho -> A rho B`

        The coefficients of these terms are typically polynomials of the
        SPAMVec's parameters, where the polynomial's variable indices index the
        *global* parameters of the SPAMVec's parent (usually a :class:`Model`)
        , not the SPAMVec's local parameter array (i.e. that returned from
        `to_vector`).


        Parameters
        ----------
        order : int
            The order of terms to get.

        return_coeff_polys : bool
            Whether a parallel list of locally-indexed (using variable indices
            corresponding to *this* object's parameters rather than its parent's)
            polynomial coefficients should be returned as well.

        Returns
        -------
        terms : list
            A list of :class:`RankOneTerm` objects.

        coefficients : list
            Only present when `return_coeff_polys == True`.
            A list of *compact* polynomial objects, meaning that each element
            is a `(vtape,ctape)` 2-tuple formed by concatenating together the
            output of :method:`Polynomial.compact`.
        """
        raise NotImplementedError("get_taylor_order_terms(...) not implemented for %s objects!" %
                                  self.__class__.__name__)

    def get_highmagnitude_terms(self, min_term_mag, force_firstorder=True, max_taylor_order=3, max_poly_vars=100):
        """
        Get the terms (from a Taylor expansion of this SPAM vector) that have
        magnitude above `min_term_mag` (the magnitude of a term is taken to
        be the absolute value of its coefficient), considering only those
        terms up to some maximum Taylor expansion order, `max_taylor_order`.

        Note that this function also *sets* the magnitudes of the returned
        terms (by calling `term.set_magnitude(...)`) based on the current
        values of this SPAM vector's parameters.  This is an essential step
        to using these terms in pruned-path-integral calculations later on.

        Parameters
        ----------
        min_term_mag : float
            the threshold for term magnitudes: only terms with magnitudes above
            this value are returned.

        force_firstorder : bool, optional
            if True, then always return all the first-order Taylor-series terms,
            even if they have magnitudes smaller than `min_term_mag`.  This
            behavior is needed for using GST with pruned-term calculations, as
            we may begin with a guess model that has no error (all terms have
            zero magnitude!) and still need to compute a meaningful jacobian at
            this point.

        max_taylor_order : int, optional
            the maximum Taylor-order to consider when checking whether term-
            magnitudes exceed `min_term_mag`.

        Returns
        -------
        highmag_terms : list
            A list of the high-magnitude terms that were found.  These
            terms are *sorted* in descending order by term-magnitude.
        first_order_indices : list
            A list of the indices into `highmag_terms` that mark which
            of these terms are first-order Taylor terms (useful when
            we're forcing these terms to always be present).
        """
        #NOTE: SAME as for LinearOperator class -- TODO consolidate in FUTURE
        #print("DB: SPAM get_high_magnitude_terms")
        v = self.to_vector()
        taylor_order = 0
        terms = []; last_len = -1; first_order_magmax = 1.0
        while len(terms) > last_len:  # while we keep adding something
            if taylor_order > 1 and first_order_magmax**taylor_order < min_term_mag:
                break  # there's no way any terms at this order reach min_term_mag - exit now!

            MAX_CACHED_TERM_ORDER = 1
            if taylor_order <= MAX_CACHED_TERM_ORDER:
                #print("order ",taylor_order," : ",len(terms), "terms")
                terms_at_order, cpolys = self.get_taylor_order_terms(taylor_order, max_poly_vars, True)
                coeffs = _bulk_eval_compact_polys_complex(
                    cpolys[0], cpolys[1], v, (len(terms_at_order),))  # an array of coeffs
                mags = _np.abs(coeffs)
                last_len = len(terms)
                #OLD: terms_at_order = [ t.copy_with_magnitude(abs(coeff)) for coeff, t in zip(coeffs, terms_at_order) ]

                if taylor_order == 1:
                    #OLD: first_order_magmax = max([t.magnitude for t in terms_at_order])
                    first_order_magmax = max(mags)

                    if force_firstorder:
                        terms.extend([(taylor_order, t.copy_with_magnitude(mag))
                                      for coeff, mag, t in zip(coeffs, mags, terms_at_order)])
                    else:
                        for mag, t in zip(mags, terms_at_order):
                            if mag >= min_term_mag:
                                terms.append((taylor_order, t.copy_with_magnitude(mag)))
                else:
                    for mag, t in zip(mags, terms_at_order):
                        if mag >= min_term_mag:
                            terms.append((taylor_order, t.copy_with_magnitude(mag)))

            else:
                terms.extend([(taylor_order, t) for t in
                              self.get_taylor_order_terms_above_mag(taylor_order,
                                                                    max_poly_vars, min_term_mag)])

            taylor_order += 1
            if taylor_order > max_taylor_order: break

        #Sort terms based on magnitude
        sorted_terms = sorted(terms, key=lambda t: t[1].magnitude, reverse=True)
        first_order_indices = [i for i, t in enumerate(sorted_terms) if t[0] == 1]
        return [t[1] for t in sorted_terms], first_order_indices

    def get_taylor_order_terms_above_mag(self, order, max_poly_vars, min_term_mag):
        """ TODO: docstring """
        v = self.to_vector()
        terms_at_order, cpolys = self.get_taylor_order_terms(order, max_poly_vars, True)
        coeffs = _bulk_eval_compact_polys_complex(
            cpolys[0], cpolys[1], v, (len(terms_at_order),))  # an array of coeffs
        terms_at_order = [t.copy_with_magnitude(abs(coeff)) for coeff, t in zip(coeffs, terms_at_order)]
        return [t for t in terms_at_order if t.magnitude >= min_term_mag]

    def frobeniusdist2(self, otherSpamVec, typ, transform=None,
                       inv_transform=None):
        """
        Return the squared frobenius difference between this spam vector and
        `otherSpamVec`, optionally transforming this vector first using
        `transform` and `inv_transform` (depending on the value of `typ`).

        Parameters
        ----------
        otherSpamVec : SPAMVec
            The other spam vector

        typ : { 'prep', 'effect' }
            Which type of SPAM vector is being transformed.

        transform, inv_transform : numpy.ndarray
            The transformation (if not None) to be performed.

        Returns
        -------
        float
        """
        vec = self.todense()
        if typ == 'prep':
            if inv_transform is None:
                return _gt.frobeniusdist2(vec, otherSpamVec.todense())
            else:
                return _gt.frobeniusdist2(_np.dot(inv_transform, vec),
                                          otherSpamVec.todense())
        elif typ == "effect":
            if transform is None:
                return _gt.frobeniusdist2(vec, otherSpamVec.todense())
            else:
                return _gt.frobeniusdist2(_np.dot(_np.transpose(transform),
                                                  vec), otherSpamVec.todense())
        else: raise ValueError("Invalid 'typ' argument: %s" % typ)

    def residuals(self, otherSpamVec, typ, transform=None, inv_transform=None):
        """
        Return a vector of residuals between this spam vector and
        `otherSpamVec`, optionally transforming this vector first using
        `transform` and `inv_transform` (depending on the value of `typ`).

        Parameters
        ----------
        otherSpamVec : SPAMVec
            The other spam vector

        typ : { 'prep', 'effect' }
            Which type of SPAM vector is being transformed.

        transform, inv_transform : numpy.ndarray
            The transformation (if not None) to be performed.

        Returns
        -------
        float
        """
        vec = self.todense()
        if typ == 'prep':
            if inv_transform is None:
                return _gt.residuals(vec, otherSpamVec.todense())
            else:
                return _gt.residuals(_np.dot(inv_transform, vec),
                                     otherSpamVec.todense())
        elif typ == "effect":
            if transform is None:
                return _gt.residuals(vec, otherSpamVec.todense())
            else:
                return _gt.residuals(_np.dot(_np.transpose(transform),
                                             vec), otherSpamVec.todense())

    def transform(self, S, typ):
        """
        Update SPAM (column) vector V as inv(S) * V or S^T * V for preparation
        or  effect SPAM vectors, respectively.

        Note that this is equivalent to state preparation vectors getting
        mapped: `rho -> inv(S) * rho` and the *transpose* of effect vectors
        being mapped as `E^T -> E^T * S`.

        Generally, the transform function updates the *parameters* of
        the SPAM vector such that the resulting vector is altered as
        described above.  If such an update cannot be done (because
        the gate parameters do not allow for it), ValueError is raised.

        Parameters
        ----------
        S : GaugeGroupElement
            A gauge group element which specifies the "S" matrix
            (and it's inverse) used in the above similarity transform.

        typ : { 'prep', 'effect' }
            Which type of SPAM vector is being transformed (see above).
        """
        if typ == 'prep':
            Si = S.get_transform_matrix_inverse()
            self.set_value(_np.dot(Si, self.todense()))
        elif typ == 'effect':
            Smx = S.get_transform_matrix()
            self.set_value(_np.dot(_np.transpose(Smx), self.todense()))
            #Evec^T --> ( Evec^T * S )^T
        else:
            raise ValueError("Invalid typ argument: %s" % typ)

    def depolarize(self, amount):
        """
        Depolarize this SPAM vector by the given `amount`.

        Generally, the depolarize function updates the *parameters* of
        the SPAMVec such that the resulting vector is depolarized.  If
        such an update cannot be done (because the gate parameters do not
        allow for it), ValueError is raised.

        Parameters
        ----------
        amount : float or tuple
            The amount to depolarize by.  If a tuple, it must have length
            equal to one less than the dimension of the gate. All but the
            first element of the spam vector (often corresponding to the
            identity element) are multiplied by `amount` (if a float) or
            the corresponding `amount[i]` (if a tuple).

        Returns
        -------
        None
        """
        if isinstance(amount, float) or _compat.isint(amount):
            D = _np.diag([1] + [1 - amount] * (self.dim - 1))
        else:
            assert(len(amount) == self.dim - 1)
            D = _np.diag([1] + list(1.0 - _np.array(amount, 'd')))
        self.set_value(_np.dot(D, self.todense()))

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return 0  # no parameters

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        return _np.array([], 'd')  # no parameters

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        assert(len(v) == 0)  # should be no parameters, and nothing to do

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.
        An empty 2D array in the StaticSPAMVec case (num_params == 0).

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        dtype = complex if self._evotype == 'statevec' else 'd'
        derivMx = _np.zeros((self.dim, 0), dtype)
        if wrtFilter is None:
            return derivMx
        else:
            return _np.take(derivMx, wrtFilter, axis=1)

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        #Default: assume Hessian can be nonzero if there are any parameters
        return self.num_params() > 0

    def hessian_wrt_params(self, wrtFilter1=None, wrtFilter2=None):
        """
        Construct the Hessian of this SPAM vector with respect to its parameters.

        This function returns a tensor whose first axis corresponds to the
        flattened operation matrix and whose 2nd and 3rd axes correspond to the
        parameters that are differentiated with respect to.

        Parameters
        ----------
        wrtFilter1, wrtFilter2 : list
            Lists of indices of the paramters to take first and second
            derivatives with respect to.  If None, then derivatives are
            taken with respect to all of the vectors's parameters.

        Returns
        -------
        numpy array
            Hessian with shape (dimension, num_params1, num_params2)
        """
        if not self.has_nonzero_hessian():
            return _np.zeros(self.size, self.num_params(), self.num_params())

        # FUTURE: create a finite differencing hessian method?
        raise NotImplementedError("hessian_wrt_params(...) is not implemented for %s objects" % self.__class__.__name__)

    #Note: no __str__ fn

    @staticmethod
    def convert_to_vector(V):
        """
        Static method that converts a vector-like object to a 2D numpy
        dim x 1 column array.

        Parameters
        ----------
        V : array_like

        Returns
        -------
        numpy array
        """
        if isinstance(V, SPAMVec):
            vector = V.todense().copy()
            vector.shape = (vector.size, 1)
        elif isinstance(V, _np.ndarray):
            vector = V.copy()
            if len(vector.shape) == 1:  # convert (N,) shape vecs to (N,1)
                vector.shape = (vector.size, 1)
        else:
            try:
                len(V)
                # XXX this is an abuse of exception handling
            except:
                raise ValueError("%s doesn't look like an array/list" % V)
            try:
                d2s = [len(row) for row in V]
            except TypeError:  # thrown if len(row) fails because no 2nd dim
                d2s = None

            if d2s is not None:
                if any([len(row) != 1 for row in V]):
                    raise ValueError("%s is 2-dimensional but 2nd dim != 1" % V)

                typ = 'd' if _np.all(_np.isreal(V)) else 'complex'
                try:
                    vector = _np.array(V, typ)  # vec is already a 2-D column vector
                except TypeError:
                    raise ValueError("%s doesn't look like an array/list" % V)
            else:
                typ = 'd' if _np.all(_np.isreal(V)) else 'complex'
                vector = _np.array(V, typ)[:, None]  # make into a 2-D column vec

        assert(len(vector.shape) == 2 and vector.shape[1] == 1)
        return vector.flatten()  # HACK for convention change -> (N,) instead of (N,1)


class DenseSPAMVec(SPAMVec):
    """
    Excapulates a parameterization of a state preparation OR POVM effect
    vector. This class is the  common base class for all specific
    parameterizations of a SPAM vector.
    """

    def __init__(self, vec, evotype, prep_or_effect):
        """ Initialize a new SPAM Vector """
        dtype = complex if evotype == "statevec" else 'd'
        vec = _np.asarray(vec, dtype=dtype)
        vec.shape = (vec.size,)  # just store 1D array flatten
        vec = _np.require(vec, requirements=['OWNDATA', 'C_CONTIGUOUS'])

        if prep_or_effect == "prep":
            if evotype == "statevec":
                rep = replib.SVStateRep(vec)
            elif evotype == "densitymx":
                rep = replib.DMStateRep(vec)
            else:
                raise ValueError("Invalid evotype for DenseSPAMVec: %s" % evotype)
        elif prep_or_effect == "effect":
            if evotype == "statevec":
                rep = replib.SVEffectRep_Dense(vec)
            elif evotype == "densitymx":
                rep = replib.DMEffectRep_Dense(vec)
            else:
                raise ValueError("Invalid evotype for DenseSPAMVec: %s" % evotype)
        else:
            raise ValueError("Invalid `prep_or_effect` argument: %s" % prep_or_effect)

        super(DenseSPAMVec, self).__init__(rep, evotype, prep_or_effect)
        assert(self.base1D.flags['C_CONTIGUOUS'] and self.base1D.flags['OWNDATA'])

    def todense(self, scratch=None):
        """
        Return this SPAM vector as a (dense) numpy array.  The memory
        in `scratch` maybe used when it is not-None.
        """
        #don't use scratch since we already have memory allocated
        return self.base1D  # *must* be a numpy array for Cython arg conversion

    @property
    def base1D(self):
        return self._rep.base

    @property
    def base(self):
        bv = self.base1D.view()
        bv.shape = (bv.size, 1)  # 'base' is by convention a (N,1)-shaped array
        return bv

    def __copy__(self):
        # We need to implement __copy__ because we defer all non-existing
        # attributes to self.base (a numpy array) which *has* a __copy__
        # implementation that we don't want to use, as it results in just a
        # copy of the numpy array.
        cls = self.__class__
        cpy = cls.__new__(cls)
        cpy.__dict__.update(self.__dict__)
        return cpy

    def __deepcopy__(self, memo):
        # We need to implement __deepcopy__ because we defer all non-existing
        # attributes to self.base (a numpy array) which *has* a __deepcopy__
        # implementation that we don't want to use, as it results in just a
        # copy of the numpy array.
        cls = self.__class__
        cpy = cls.__new__(cls)
        memo[id(self)] = cpy
        for k, v in self.__dict__.items():
            setattr(cpy, k, _copy.deepcopy(v, memo))
        return cpy

    #Access to underlying array
    def __getitem__(self, key):
        self.dirty = True
        return self.base.__getitem__(key)

    def __getslice__(self, i, j):
        self.dirty = True
        return self.__getitem__(slice(i, j))  # Called for A[:]

    def __setitem__(self, key, val):
        self.dirty = True
        return self.base.__setitem__(key, val)

    def __getattr__(self, attr):
        #use __dict__ so no chance for recursive __getattr__
        if '_rep' in self.__dict__:  # sometimes in loading __getattr__ gets called before the instance is loaded
            ret = getattr(self.base, attr)
        else:
            raise AttributeError("No attribute:", attr)
        self.dirty = True
        return ret

    #Mimic array
    def __pos__(self): return self.base
    def __neg__(self): return -self.base
    def __abs__(self): return abs(self.base)
    def __add__(self, x): return self.base + x
    def __radd__(self, x): return x + self.base
    def __sub__(self, x): return self.base - x
    def __rsub__(self, x): return x - self.base
    def __mul__(self, x): return self.base * x
    def __rmul__(self, x): return x * self.base
    def __truediv__(self, x): return self.base / x
    def __rtruediv__(self, x): return x / self.base
    def __floordiv__(self, x): return self.base // x
    def __rfloordiv__(self, x): return x // self.base
    def __pow__(self, x): return self.base ** x
    def __eq__(self, x): return self.base == x
    def __len__(self): return len(self.base)
    def __int__(self): return int(self.base)
    def __long__(self): return int(self.base)
    def __float__(self): return float(self.base)
    def __complex__(self): return complex(self.base)

    def __str__(self):
        s = "%s with dimension %d\n" % (self.__class__.__name__, self.dim)
        s += _mt.mx_to_string(self.todense(), width=4, prec=2)
        return s


class StaticSPAMVec(DenseSPAMVec):
    """
    Encapsulates a SPAM vector that is completely fixed, or "static", meaning
      that is contains no parameters.
    """

    def __init__(self, vec, evotype="auto", typ="prep"):
        """
        Initialize a StaticSPAMVec object.

        Parameters
        ----------
        vec : array_like or SPAMVec
            a 1D numpy array representing the SPAM operation.  The
            shape of this array sets the dimension of the SPAM op.

        evotype : {"densitymx", "statevec"}
            the evolution type being used.

        typ : {"prep", "effect"}
            whether this is a state preparation or an effect (measurement)
            SPAM vector.
        """
        vec = SPAMVec.convert_to_vector(vec)
        if evotype == "auto":
            evotype = "statevec" if _np.iscomplexobj(vec) else "densitymx"
        elif evotype == "statevec":
            vec = _np.asarray(vec, complex)  # ensure all statevec vecs are complex (densitymx could be either?)

        assert(evotype in ("statevec", "densitymx")), \
            "Invalid evolution type '%s' for %s" % (evotype, self.__class__.__name__)
        DenseSPAMVec.__init__(self, vec, evotype, typ)


class FullSPAMVec(DenseSPAMVec):
    """
    Encapsulates a SPAM vector that is fully parameterized, that is,
      each element of the SPAM vector is an independent parameter.
    """

    def __init__(self, vec, evotype="auto", typ="prep"):
        """
        Initialize a FullSPAMVec object.

        Parameters
        ----------
        vec : array_like or SPAMVec
            a 1D numpy array representing the SPAM operation.  The
            shape of this array sets the dimension of the SPAM op.

        evotype : {"densitymx", "statevec"}
            the evolution type being used.

        typ : {"prep", "effect"}
            whether this is a state preparation or an effect (measurement)
            SPAM vector.
        """
        vec = SPAMVec.convert_to_vector(vec)
        if evotype == "auto":
            evotype = "statevec" if _np.iscomplexobj(vec) else "densitymx"
        assert(evotype in ("statevec", "densitymx")), \
            "Invalid evolution type '%s' for %s" % (evotype, self.__class__.__name__)
        DenseSPAMVec.__init__(self, vec, evotype, typ)

    def set_value(self, vec):
        """
        Attempts to modify SPAMVec parameters so that the specified raw
        SPAM vector becomes vec.  Will raise ValueError if this operation
        is not possible.

        Parameters
        ----------
        vec : array_like or SPAMVec
            A numpy array representing a SPAM vector, or a SPAMVec object.

        Returns
        -------
        None
        """
        vec = SPAMVec.convert_to_vector(vec)
        if(vec.size != self.dim):
            raise ValueError("Argument must be length %d" % self.dim)
        self.base1D[:] = vec
        self.dirty = True

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return 2 * self.size if self._evotype == "statevec" else self.size

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        if self._evotype == "statevec":
            return _np.concatenate((self.base1D.real, self.base1D.imag), axis=0)
        else:
            return self.base1D

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        if self._evotype == "statevec":
            self.base1D[:] = v[0:self.dim] + 1j * v[self.dim:]
        else:
            self.base1D[:] = v
        if not nodirty: self.dirty = True

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        if self._evotype == "statevec":
            derivMx = _np.concatenate((_np.identity(self.dim, complex),
                                       1j * _np.identity(self.dim, complex)), axis=1)
        else:
            derivMx = _np.identity(self.dim, 'd')

        if wrtFilter is None:
            return derivMx
        else:
            return _np.take(derivMx, wrtFilter, axis=1)

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        return False


class TPSPAMVec(DenseSPAMVec):
    """
    Encapsulates a SPAM vector that is fully parameterized except for the first
    element, which is frozen to be 1/(d**0.25).  This is so that, when the SPAM
    vector is interpreted in the Pauli or Gell-Mann basis, the represented
    density matrix has trace == 1.  This restriction is frequently used in
    conjuction with trace-preserving (TP) gates.
    """

    #Note: here we assume that the first basis element is (1/sqrt(x) * I),
    # where I the d-dimensional identity (where len(vector) == d**2). So
    # if Tr(basisEl*basisEl) == Tr(1/x*I) == d/x must == 1, then we must
    # have x == d.  Thus, we multiply this first basis element by
    # alpha = 1/sqrt(d) to obtain a trace-1 matrix, i.e., finding alpha
    # s.t. Tr(alpha*[1/sqrt(d)*I]) == 1 => alpha*d/sqrt(d) == 1 =>
    # alpha = 1/sqrt(d) = 1/(len(vec)**0.25).
    def __init__(self, vec):
        """
        Initialize a TPSPAMVec object.

        Parameters
        ----------
        vec : array_like or SPAMVec
            a 1D numpy array representing the SPAM operation.  The
            shape of this array sets the dimension of the SPAM op.
        """
        vector = SPAMVec.convert_to_vector(vec)
        firstEl = len(vector)**-0.25
        if not _np.isclose(vector[0], firstEl):
            raise ValueError("Cannot create TPSPAMVec: "
                             "first element must equal %g!" % firstEl)
        DenseSPAMVec.__init__(self, vec, "densitymx", "prep")
        assert(isinstance(self.base, _ProtectedArray))

    @property
    def base(self):
        bv = self.base1D.view()
        bv.shape = (bv.size, 1)
        return _ProtectedArray(bv, indicesToProtect=(0, 0))

    def set_value(self, vec):
        """
        Attempts to modify SPAMVec parameters so that the specified raw
        SPAM vector becomes vec.  Will raise ValueError if this operation
        is not possible.

        Parameters
        ----------
        vec : array_like or SPAMVec
            A numpy array representing a SPAM vector, or a SPAMVec object.

        Returns
        -------
        None
        """
        vec = SPAMVec.convert_to_vector(vec)
        firstEl = (self.dim)**-0.25
        if(vec.size != self.dim):
            raise ValueError("Argument must be length %d" % self.dim)
        if not _np.isclose(vec[0], firstEl):
            raise ValueError("Cannot create TPSPAMVec: "
                             "first element must equal %g!" % firstEl)
        self.base1D[1:] = vec[1:]
        self.dirty = True

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return self.dim - 1

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        return self.base1D[1:]  # .real in case of complex matrices?

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        assert(_np.isclose(self.base1D[0], (self.dim)**-0.25))
        self.base1D[1:] = v
        if not nodirty: self.dirty = True

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        derivMx = _np.identity(self.dim, 'd')  # TP vecs assumed real
        derivMx = derivMx[:, 1:]  # remove first col ( <=> first-el parameters )
        if wrtFilter is None:
            return derivMx
        else:
            return _np.take(derivMx, wrtFilter, axis=1)

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        return False


class ComplementSPAMVec(DenseSPAMVec):
    """
    Encapsulates a SPAM vector that is parameterized by
    `I - sum(other_spam_vecs)` where `I` is a (static) identity element
    and `other_param_vecs` is a list of other spam vectors in the same parent
    :class:`POVM`.  This only *partially* implements the SPAMVec interface
    (some methods such as `to_vector` and `from_vector` will thunk down to base
    class versions which raise `NotImplementedError`), as instances are meant to
    be contained within a :class:`POVM` which takes care of vectorization.
    """

    def __init__(self, identity, other_spamvecs):
        """
        Initialize a ComplementSPAMVec object.

        Parameters
        ----------
        identity : array_like or SPAMVec
            a 1D numpy array representing the static identity operation from
            which the sum of the other vectors is subtracted.

        other_spamvecs : list of SPAMVecs
            A list of the "other" parameterized SPAM vectors which are
            subtracted from `identity` to compute the final value of this
            "complement" SPAM vector.
        """
        self.identity = FullSPAMVec(
            SPAMVec.convert_to_vector(identity))  # so easy to transform
        # or depolarize by parent POVM

        self.other_vecs = other_spamvecs
        #Note: we assume that our parent will do the following:
        # 1) set our gpindices to indicate how many parameters we have
        # 2) set the gpindices of the elements of other_spamvecs so
        #    that they index into our local parameter vector.

        DenseSPAMVec.__init__(self, self.identity, "densitymx", "effect")  # dummy
        self._construct_vector()  # reset's self.base

    def _construct_vector(self):
        self.base1D.flags.writeable = True
        self.base1D[:] = self.identity.base1D - sum([vec.base1D for vec in self.other_vecs])
        self.base1D.flags.writeable = False

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return len(self.gpindices_as_array())

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        raise ValueError(("ComplementSPAMVec.to_vector() should never be called"
                          " - use TPPOVM.to_vector() instead"))

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize this SPAM vector using a vector of its parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of parameters.

        Returns
        -------
        None
        """
        #Rely on prior .from_vector initialization of self.other_vecs, so
        # we just construct our vector based on them.
        #Note: this is needed for finite-differencing in map-based calculator
        self._construct_vector()
        if not nodirty: self.dirty = True

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        if len(self.other_vecs) == 0: return _np.zeros((self.dim, 0), 'd')  # Complement vecs assumed real
        Np = len(self.gpindices_as_array())
        neg_deriv = _np.zeros((self.dim, Np), 'd')
        for ovec in self.other_vecs:
            local_inds = _modelmember._decompose_gpindices(
                self.gpindices, ovec.gpindices)
            #Note: other_vecs are not copies but other *sibling* effect vecs
            # so their gpindices index the same space as this complement vec's
            # does - so we need to "_decompose_gpindices"
            neg_deriv[:, local_inds] += ovec.deriv_wrt_params()
        derivMx = -neg_deriv

        if wrtFilter is None:
            return derivMx
        else:
            return _np.take(derivMx, wrtFilter, axis=1)

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        return False


class CPTPSPAMVec(DenseSPAMVec):
    """
    Encapsulates a SPAM vector that is parameterized through the Cholesky
    decomposition of it's standard-basis representation as a density matrix
    (not a Liouville vector).  The resulting SPAM vector thus represents a
    positive density matrix, and additional constraints on the parameters
    also guarantee that the trace == 1.  This SPAM vector is meant for
    use with CPTP processes, hence the name.
    """

    def __init__(self, vec, basis, truncate=False):
        """
        Initialize a CPTPSPAMVec object.

        Parameters
        ----------
        vec : array_like or SPAMVec
            a 1D numpy array representing the SPAM operation.  The
            shape of this array sets the dimension of the SPAM op.

        basis : {"std", "gm", "pp", "qt"} or Basis
            The basis `vec` is in.  Needed because this parameterization
            requires we construct the density matrix corresponding to
            the Lioville vector `vec`.

        trunctate : bool, optional
            Whether or not a non-positive, trace=1 `vec` should
            be truncated to force a successful construction.
        """
        vector = SPAMVec.convert_to_vector(vec)
        basis = _Basis.cast(basis, len(vector))

        self.basis = basis
        self.basis_mxs = basis.elements  # shape (len(vec), dmDim, dmDim)
        self.basis_mxs = _np.rollaxis(self.basis_mxs, 0, 3)  # shape (dmDim, dmDim, len(vec))
        assert(self.basis_mxs.shape[-1] == len(vector))

        # set self.params and self.dmDim
        self._set_params_from_vector(vector, truncate)

        #scratch space
        self.Lmx = _np.zeros((self.dmDim, self.dmDim), 'complex')

        DenseSPAMVec.__init__(self, vector, "densitymx", "prep")

    def _set_params_from_vector(self, vector, truncate):
        density_mx = _np.dot(self.basis_mxs, vector)
        density_mx = density_mx.squeeze()
        dmDim = density_mx.shape[0]
        assert(dmDim == density_mx.shape[1]), "Density matrix must be square!"

        trc = _np.trace(density_mx)
        assert(truncate or _np.isclose(trc, 1.0)), \
            "`vec` must correspond to a trace-1 density matrix (truncate == False)!"

        if not _np.isclose(trc, 1.0):  # truncate to trace == 1
            density_mx -= _np.identity(dmDim, 'd') / dmDim * (trc - 1.0)

        #push any slightly negative evals of density_mx positive
        # so that the Cholesky decomp will work.
        evals, U = _np.linalg.eig(density_mx)
        Ui = _np.linalg.inv(U)

        assert(truncate or all([ev >= -1e-12 for ev in evals])), \
            "`vec` must correspond to a positive density matrix (truncate == False)!"

        pos_evals = evals.clip(1e-16, 1e100)
        density_mx = _np.dot(U, _np.dot(_np.diag(pos_evals), Ui))
        try:
            Lmx = _np.linalg.cholesky(density_mx)
        except _np.linalg.LinAlgError:  # Lmx not postitive definite?
            pos_evals = evals.clip(1e-12, 1e100)  # try again with 1e-12
            density_mx = _np.dot(U, _np.dot(_np.diag(pos_evals), Ui))
            Lmx = _np.linalg.cholesky(density_mx)

        #check TP condition: that diagonal els of Lmx squared add to 1.0
        Lmx_norm = _np.trace(_np.dot(Lmx.T.conjugate(), Lmx))  # sum of magnitude^2 of all els
        assert(_np.isclose(Lmx_norm, 1.0)), \
            "Cholesky decomp didn't preserve trace=1!"

        self.dmDim = dmDim
        self.params = _np.empty(dmDim**2, 'd')
        for i in range(dmDim):
            assert(_np.linalg.norm(_np.imag(Lmx[i, i])) < IMAG_TOL)
            self.params[i * dmDim + i] = Lmx[i, i].real  # / paramNorm == 1 as asserted above
            for j in range(i):
                self.params[i * dmDim + j] = Lmx[i, j].real
                self.params[j * dmDim + i] = Lmx[i, j].imag

    def _construct_vector(self):
        dmDim = self.dmDim

        #  params is an array of length dmDim^2-1 that
        #  encodes a lower-triangular matrix "Lmx" via:
        #  Lmx[i,i] = params[i*dmDim + i] / param-norm  # i = 0...dmDim-2
        #     *last diagonal el is given by sqrt(1.0 - sum(L[i,j]**2))
        #  Lmx[i,j] = params[i*dmDim + j] + 1j*params[j*dmDim+i] (i > j)

        param2Sum = _np.vdot(self.params, self.params)  # or "dot" would work, since params are real
        paramNorm = _np.sqrt(param2Sum)  # also the norm of *all* Lmx els

        for i in range(dmDim):
            self.Lmx[i, i] = self.params[i * dmDim + i] / paramNorm
            for j in range(i):
                self.Lmx[i, j] = (self.params[i * dmDim + j] + 1j * self.params[j * dmDim + i]) / paramNorm

        Lmx_norm = _np.trace(_np.dot(self.Lmx.T.conjugate(), self.Lmx))  # sum of magnitude^2 of all els
        assert(_np.isclose(Lmx_norm, 1.0)), "Violated trace=1 condition!"

        #The (complex, Hermitian) density matrix is build by
        # assuming Lmx is its Cholesky decomp, which makes
        # the density matrix is pos-def.
        density_mx = _np.dot(self.Lmx, self.Lmx.T.conjugate())
        assert(_np.isclose(_np.trace(density_mx), 1.0)), "density matrix must be trace == 1"

        # write density matrix in given basis: = sum_i alpha_i B_i
        # ASSUME that basis is orthogonal, i.e. Tr(Bi^dag*Bj) = delta_ij
        basis_mxs = _np.rollaxis(self.basis_mxs, 2)  # shape (dmDim, dmDim, len(vec))
        vec = _np.array([_np.trace(_np.dot(M.T.conjugate(), density_mx)) for M in basis_mxs])

        #for now, assume Liouville vector should always be real (TODO: add 'real' flag later?)
        assert(_np.linalg.norm(_np.imag(vec)) < IMAG_TOL)
        vec = _np.real(vec)

        self.base1D.flags.writeable = True
        self.base1D[:] = vec[:]  # so shape is (dim,1) - the convention for spam vectors
        self.base1D.flags.writeable = False

    def set_value(self, vec):
        """
        Attempts to modify SPAMVec parameters so that the specified raw
        SPAM vector becomes vec.  Will raise ValueError if this operation
        is not possible.

        Parameters
        ----------
        vec : array_like or SPAMVec
            A numpy array representing a SPAM vector, or a SPAMVec object.

        Returns
        -------
        None
        """
        try:
            self._set_params_from_vector(vec, truncate=False)
            self.dirty = True
        except AssertionError as e:
            raise ValueError("Error initializing the parameters of this "
                             "CPTPSPAMVec object: " + str(e))

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        assert(self.dmDim**2 == self.dim)  # should at least be true without composite bases...
        return self.dmDim**2

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        return self.params

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        assert(len(v) == self.num_params())
        self.params[:] = v[:]
        self._construct_vector()
        if not nodirty: self.dirty = True

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        dmDim = self.dmDim
        nP = len(self.params)
        assert(nP == dmDim**2)  # number of parameters

        # v_i = trace( B_i^dag * Lmx * Lmx^dag )
        # d(v_i) = trace( B_i^dag * (dLmx * Lmx^dag + Lmx * (dLmx)^dag) )  #trace = linear so commutes w/deriv
        #               /
        # where dLmx/d[ab] = {
        #               \
        L, Lbar = self.Lmx, self.Lmx.conjugate()
        F1 = _np.tril(_np.ones((dmDim, dmDim), 'd'))
        F2 = _np.triu(_np.ones((dmDim, dmDim), 'd'), 1) * 1j
        conj_basis_mxs = self.basis_mxs.conjugate()

        # Derivative of vector wrt params; shape == [vecLen,dmDim,dmDim] *not dealing with TP condition yet*
        # (first get derivative assuming last diagonal el of Lmx *is* a parameter, then use chain rule)
        dVdp = _np.einsum('aml,mb,ab->lab', conj_basis_mxs, Lbar, F1)  # only a >= b nonzero (F1)
        dVdp += _np.einsum('mal,mb,ab->lab', conj_basis_mxs, L, F1)    # ditto
        dVdp += _np.einsum('bml,ma,ab->lab', conj_basis_mxs, Lbar, F2)  # only b > a nonzero (F2)
        dVdp += _np.einsum('mbl,ma,ab->lab', conj_basis_mxs, L, F2.conjugate())  # ditto

        dVdp.shape = [dVdp.shape[0], nP]  # jacobian with respect to "p" params,
        # which don't include normalization for TP-constraint

        #Now get jacobian of actual params wrt the params used above. Denote the actual
        # params "P" in variable names, so p_ij = P_ij / sqrt(sum(P_xy**2))
        param2Sum = _np.vdot(self.params, self.params)
        paramNorm = _np.sqrt(param2Sum)  # norm of *all* Lmx els (note lastDiagEl
        dpdP = _np.identity(nP, 'd')

        # all p_ij params ==  P_ij / paramNorm = P_ij / sqrt(sum(P_xy**2))
        # and so have derivs wrt *all* Pxy elements.
        for ij in range(nP):
            for kl in range(nP):
                if ij == kl:
                    # dp_ij / dP_ij = 1.0 / (sum(P_xy**2))^(1/2) - 0.5 * P_ij / (sum(P_xy**2))^(3/2) * 2*P_ij
                    #               = 1.0 / (sum(P_xy**2))^(1/2) - P_ij^2 / (sum(P_xy**2))^(3/2)
                    dpdP[ij, ij] = 1.0 / paramNorm - self.params[ij]**2 / paramNorm**3
                else:
                    # dp_ij / dP_kl = -0.5 * P_ij / (sum(P_xy**2))^(3/2) * 2*P_kl
                    #               = - P_ij * P_kl / (sum(P_xy**2))^(3/2)
                    dpdP[ij, kl] = - self.params[ij] * self.params[kl] / paramNorm**3

        #Apply the chain rule to get dVdP:
        dVdP = _np.dot(dVdp, dpdP)  # shape (vecLen, nP) - the jacobian!
        dVdp = dpdP = None  # free memory!

        assert(_np.linalg.norm(_np.imag(dVdP)) < IMAG_TOL)
        derivMx = _np.real(dVdP)

        if wrtFilter is None:
            return derivMx
        else:
            return _np.take(derivMx, wrtFilter, axis=1)

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        return True

    def hessian_wrt_params(self, wrtFilter1=None, wrtFilter2=None):
        """
        Construct the Hessian of this SPAM vector with respect to its parameters.

        This function returns a tensor whose first axis corresponds to the
        flattened operation matrix and whose 2nd and 3rd axes correspond to the
        parameters that are differentiated with respect to.

        Parameters
        ----------
        wrtFilter1, wrtFilter2 : list
            Lists of indices of the paramters to take first and second
            derivatives with respect to.  If None, then derivatives are
            taken with respect to all of the vectors's parameters.

        Returns
        -------
        numpy array
            Hessian with shape (dimension, num_params1, num_params2)
        """
        raise NotImplementedError("TODO: add hessian computation for CPTPSPAMVec")


class TensorProdSPAMVec(SPAMVec):
    """
    Encapsulates a SPAM vector that is a tensor-product of other SPAM vectors.
    """

    def __init__(self, typ, factors, povmEffectLbls=None):
        """
        Initialize a TensorProdSPAMVec object.

        Parameters
        ----------
        typ : {"prep","effect"}
            The type of spam vector - either a product of preparation
            vectors ("prep") or of POVM effect vectors ("effect")

        factors : list of SPAMVecs or POVMs
            if `typ == "prep"`, a list of the component SPAMVecs; if
            `typ == "effect"`, a list of "reference" POVMs into which
            `povmEffectLbls` indexes.

        povmEffectLbls : array-like
            Only non-None when `typ == "effect"`.  The effect label of each
            factor POVM which is tensored together to form this SPAM vector.
        """
        assert(len(factors) > 0), "Must have at least one factor!"

        self.factors = factors  # do *not* copy - needs to reference common objects
        self.Np = sum([fct.num_params() for fct in factors])
        if typ == "effect":
            self.effectLbls = _np.array(povmEffectLbls)
        elif typ == "prep":
            assert(povmEffectLbls is None), '`povmEffectLbls` must be None when `typ != "effects"`'
            self.effectLbls = None
        else: raise ValueError("Invalid `typ` argument: %s" % typ)

        dim = _np.product([fct.dim for fct in factors])
        evotype = self.factors[0]._evotype

        #Arrays for speeding up kron product in effect reps
        if evotype in ("statevec", "densitymx") and typ == "effect":  # types that require fast kronecker prods
            max_factor_dim = max(fct.dim for fct in factors)
            self._fast_kron_array = _np.ascontiguousarray(
                _np.empty((len(factors), max_factor_dim), complex if evotype == "statevec" else 'd'))
            self._fast_kron_factordims = _np.ascontiguousarray(
                _np.array([fct.dim for fct in factors], _np.int64))
        else:  # "stabilizer", "svterm", "cterm"
            self._fast_kron_array = None
            self._fast_kron_factordims = None

        #Create representation
        if evotype == "statevec":
            if typ == "prep":  # prep-type vectors can be represented as dense effects too
                rep = replib.SVStateRep(_np.ascontiguousarray(_np.zeros(dim, complex)))
                #sometimes: return replib.SVEffectRep_Dense(self.todense()) ???
            else:  # "effect"
                rep = replib.SVEffectRep_TensorProd(self._fast_kron_array, self._fast_kron_factordims,
                                                    len(self.factors), self._fast_kron_array.shape[1], dim)
        elif evotype == "densitymx":
            if typ == "prep":
                vec = _np.require(_np.zeros(dim, 'd'), requirements=['OWNDATA', 'C_CONTIGUOUS'])
                rep = replib.DMStateRep(vec)
                #sometimes: return replib.DMEffectRep_Dense(self.todense()) ???
            else:  # "effect"
                rep = replib.DMEffectRep_TensorProd(self._fast_kron_array, self._fast_kron_factordims,
                                                    len(self.factors), self._fast_kron_array.shape[1], dim)
        elif evotype == "stabilizer":
            if typ == "prep":
                #Rep is stabilizer-rep tuple, just like StabilizerSPAMVec
                sframe_factors = [f.todense() for f in self.factors]  # StabilizerFrame for each factor
                rep = _stabilizer.sframe_kronecker(sframe_factors).torep()

                #Sometimes ???
                # prep-type vectors can be represented as dense effects too; this just means that self.factors
                # => self.factors should all be StabilizerSPAMVec objs
                #else:  # self._prep_or_effect == "effect", so each factor is a StabilizerEffectVec
                #  outcomes = _np.array(list(_itertools.chain(*[f.outcomes for f in self.factors])), _np.int64)
                #  return replib.SBEffectRep(outcomes)

            else:  # self._prep_or_effect == "effect", so each factor is a StabilizerZPOVM
                # like above, but get a StabilizerEffectVec from each StabilizerZPOVM factor
                factorPOVMs = self.factors
                factorVecs = [factorPOVMs[i][self.effectLbls[i]] for i in range(1, len(factorPOVMs))]
                outcomes = _np.array(list(_itertools.chain(*[f.outcomes for f in factorVecs])), _np.int64)
                rep = replib.SBEffectRep(outcomes)

                #OLD - now can remove outcomes prop?
                #raise ValueError("Cannot convert Stabilizer tensor product effect to a representation!")
                # should be using effect.outcomes property...

        else:  # self._evotype in ("svterm","cterm")
            rep = dim  # no reps for term-based evotypes

        SPAMVec.__init__(self, rep, evotype, typ)
        self._update_rep()  # initializes rep data
        #sets gpindices, so do before stuff below

        if typ == "effect":
            #Set our parent and gpindices based on those of factor-POVMs, which
            # should all be owned by a TensorProdPOVM object.
            # (for now say we depend on *all* the POVMs parameters (even though
            #  we really only depend on one element of each POVM, which may allow
            #  using just a subset of each factor POVMs indices - but this is tricky).
            self.set_gpindices(_slct.list_to_slice(
                _np.concatenate([fct.gpindices_as_array()
                                 for fct in factors], axis=0), True, False),
                               factors[0].parent)  # use parent of first factor
            # (they should all be the same)
        else:
            # don't init our own gpindices (prep case), since our parent
            # is likely to be a Model and it will init them correctly.
            #But do set the indices of self.factors, since they're now
            # considered "owned" by this product-prep-vec (different from
            # the "effect" case when the factors are shared).
            off = 0
            for fct in factors:
                assert(isinstance(fct, SPAMVec)), "Factors must be SPAMVec objects!"
                N = fct.num_params()
                fct.set_gpindices(slice(off, off + N), self); off += N
            assert(off == self.Np)

    def _fill_fast_kron(self):
        """ Fills in self._fast_kron_array based on current self.factors """
        if self._prep_or_effect == "prep":
            for i, factor_dim in enumerate(self._fast_kron_factordims):
                self._fast_kron_array[i][0:factor_dim] = self.factors[i].todense()
        else:
            factorPOVMs = self.factors
            for i, (factor_dim, Elbl) in enumerate(zip(self._fast_kron_factordims, self.effectLbls)):
                self._fast_kron_array[i][0:factor_dim] = factorPOVMs[i][Elbl].todense()

    def _update_rep(self):
        if self._evotype in ("statevec", "densitymx"):
            if self._prep_or_effect == "prep":
                self._rep.base[:] = self.todense()
            else:
                self._fill_fast_kron()  # updates effect reps
        elif self._evotype == "stabilizer":
            if self._prep_or_effect == "prep":
                #we need to update self._rep, which is a SBStateRep object.  For now, we
                # kinda punt and just create a new rep and copy its data over to the existing
                # rep instead of having some kind of update method within SBStateRep...
                # (TODO FUTURE - at least a .copy_from method?)
                sframe_factors = [f.todense() for f in self.factors]  # StabilizerFrame for each factor
                new_rep = _stabilizer.sframe_kronecker(sframe_factors).torep()
                self._rep.smatrix[:, :] = new_rep.smatrix[:, :]
                self._rep.pvectors[:, :] = new_rep.pvectors[:, :]
                self._rep.amps[:, :] = new_rep.amps[:, :]
            else:
                pass  # I think the below (e.g. 'outcomes') is not altered by any parameters
                #factorPOVMs = self.factors
                #factorVecs = [factorPOVMs[i][self.effectLbls[i]] for i in range(1, len(factorPOVMs))]
                #outcomes = _np.array(list(_itertools.chain(*[f.outcomes for f in factorVecs])), _np.int64)
                #rep = replib.SBEffectRep(outcomes)

    def todense(self):
        """
        Return this SPAM vector as a (dense) numpy array.
        """
        if self._evotype in ("statevec", "densitymx"):
            if len(self.factors) == 0: return _np.empty(0, complex if self._evotype == "statevec" else 'd')
            #NOTE: moved a fast version of todense to replib - could use that if need a fast todense call...

            if self._prep_or_effect == "prep":
                ret = self.factors[0].todense()  # factors are just other SPAMVecs
                for i in range(1, len(self.factors)):
                    ret = _np.kron(ret, self.factors[i].todense())
            else:
                factorPOVMs = self.factors
                ret = factorPOVMs[0][self.effectLbls[0]].todense()
                for i in range(1, len(factorPOVMs)):
                    ret = _np.kron(ret, factorPOVMs[i][self.effectLbls[i]].todense())

            return ret
        elif self._evotype == "stabilizer":

            if self._prep_or_effect == "prep":
                # => self.factors should all be StabilizerSPAMVec objs
                #Return stabilizer-rep tuple, just like StabilizerSPAMVec
                sframe_factors = [f.todense() for f in self.factors]
                return _stabilizer.sframe_kronecker(sframe_factors)

            else:  # self._prep_or_effect == "effect", so each factor is a StabilizerEffectVec
                raise ValueError("Cannot convert Stabilizer tensor product effect to an array!")
                # should be using effect.outcomes property...
        else:  # self._evotype in ("svterm","cterm")
            raise NotImplementedError("todense() not implemented for %s evolution type" % self._evotype)

    #def torep(self, typ, outrep=None):
    #    """
    #    Return a "representation" object for this SPAM vector.
    #
    #    Such objects are primarily used internally by pyGSTi to compute
    #    things like probabilities more efficiently.
    #
    #    Parameters
    #    ----------
    #    typ : {'prep','effect'}
    #        The type of representation (for cases when the vector type is
    #        not already defined).
    #
    #    outrep : StateRep
    #        If not None, an existing state representation appropriate to this
    #        SPAM vector that may be used instead of allocating a new one.
    #
    #    Returns
    #    -------
    #    StateRep
    #    """
    #    assert(len(self.factors) > 0), "Cannot get representation of a TensorProdSPAMVec with no factors!"
    #    assert(self._prep_or_effect in ('prep', 'effect')), "Invalid internal type: %s!" % self._prep_or_effect
    #
    #    #FUTURE: use outrep as scratch for rep constructor?
    #    if self._evotype == "statevec":
    #        if self._prep_or_effect == "prep":  # prep-type vectors can be represented as dense effects too
    #            if typ == "prep":
    #                return replib.SVStateRep(self.todense())
    #            else:
    #                return replib.SVEffectRep_Dense(self.todense())
    #        else:
    #            return replib.SVEffectRep_TensorProd(self._fast_kron_array, self._fast_kron_factordims,
    #                                                 len(self.factors), self._fast_kron_array.shape[1], self.dim)
    #    elif self._evotype == "densitymx":
    #        if self._prep_or_effect == "prep":
    #            if typ == "prep":
    #                return replib.DMStateRep(self.todense())
    #            else:
    #                return replib.DMEffectRep_Dense(self.todense())
    #
    #        else:
    #            return replib.DMEffectRep_TensorProd(self._fast_kron_array, self._fast_kron_factordims,
    #                                                 len(self.factors), self._fast_kron_array.shape[1], self.dim)
    #
    #    elif self._evotype == "stabilizer":
    #        if self._prep_or_effect == "prep":
    #            # prep-type vectors can be represented as dense effects too; this just means that self.factors
    #            if typ == "prep":
    #                # => self.factors should all be StabilizerSPAMVec objs
    #                #Return stabilizer-rep tuple, just like StabilizerSPAMVec
    #                sframe_factors = [f.todense() for f in self.factors]  # StabilizerFrame for each factor
    #                return _stabilizer.sframe_kronecker(sframe_factors).torep()
    #            else:  # self._prep_or_effect == "effect", so each factor is a StabilizerEffectVec
    #                outcomes = _np.array(list(_itertools.chain(*[f.outcomes for f in self.factors])), _np.int64)
    #                return replib.SBEffectRep(outcomes)
    #
    #        else:  # self._prep_or_effect == "effect", so each factor is a StabilizerZPOVM
    #            # like above, but get a StabilizerEffectVec from each StabilizerZPOVM factor
    #            factorPOVMs = self.factors
    #            factorVecs = [factorPOVMs[i][self.effectLbls[i]] for i in range(1, len(factorPOVMs))]
    #            outcomes = _np.array(list(_itertools.chain(*[f.outcomes for f in factorVecs])), _np.int64)
    #            return replib.SBEffectRep(outcomes)
    #
    #            #OLD - now can remove outcomes prop?
    #            #raise ValueError("Cannot convert Stabilizer tensor product effect to a representation!")
    #            # should be using effect.outcomes property...
    #
    #    else:  # self._evotype in ("svterm","cterm")
    #        raise NotImplementedError("torep() not implemented for %s evolution type" % self._evotype)

    def get_taylor_order_terms(self, order, max_poly_vars=100, return_coeff_polys=False):
        """
        Get the `order`-th order Taylor-expansion terms of this SPAM vector.

        This function either constructs or returns a cached list of the terms at
        the given order.  Each term is "rank-1", meaning that it is a state
        preparation followed by or POVM effect preceded by actions on a
        density matrix `rho` of the form:

        `rho -> A rho B`

        The coefficients of these terms are typically polynomials of the
        SPAMVec's parameters, where the polynomial's variable indices index the
        *global* parameters of the SPAMVec's parent (usually a :class:`Model`)
        , not the SPAMVec's local parameter array (i.e. that returned from
        `to_vector`).


        Parameters
        ----------
        order : int
            The order of terms to get.

        return_coeff_polys : bool
            Whether a parallel list of locally-indexed (using variable indices
            corresponding to *this* object's parameters rather than its parent's)
            polynomial coefficients should be returned as well.

        Returns
        -------
        terms : list
            A list of :class:`RankOneTerm` objects.

        coefficients : list
            Only present when `return_coeff_polys == True`.
            A list of *compact* polynomial objects, meaning that each element
            is a `(vtape,ctape)` 2-tuple formed by concatenating together the
            output of :method:`Polynomial.compact`.
        """
        from .operation import EmbeddedOp as _EmbeddedGateMap
        terms = []
        fnq = [int(round(_np.log2(f.dim))) // 2 for f in self.factors]  # num of qubits per factor
        # assumes density matrix evolution
        total_nQ = sum(fnq)  # total number of qubits

        for p in _lt.partition_into(order, len(self.factors)):
            if self._prep_or_effect == "prep":
                factor_lists = [self.factors[i].get_taylor_order_terms(pi, max_poly_vars) for i, pi in enumerate(p)]
            else:
                factorPOVMs = self.factors
                factor_lists = [factorPOVMs[i][Elbl].get_taylor_order_terms(pi, max_poly_vars)
                                for i, (pi, Elbl) in enumerate(zip(p, self.effectLbls))]

            # When possible, create COLLAPSED factor_lists so each factor has just a single
            # (SPAMVec) pre & post op, which can be formed into the new terms'
            # TensorProdSPAMVec ops.
            # - DON'T collapse stabilizer states & clifford ops - can't for POVMs
            collapsible = False  # bool(self._evotype =="svterm") # need to use reps for collapsing now... TODO?

            if collapsible:
                factor_lists = [[t.collapse_vec() for t in fterms] for fterms in factor_lists]

            for factors in _itertools.product(*factor_lists):
                # create a term with a TensorProdSPAMVec - Note we always create
                # "prep"-mode vectors, since even when self._prep_or_effect == "effect" these
                # vectors are created with factor (prep- or effect-type) SPAMVecs not factor POVMs
                # we workaround this by still allowing such "prep"-mode
                # TensorProdSPAMVecs to be represented as effects (i.e. in torep('effect'...) works)
                coeff = _functools.reduce(lambda x, y: x.mult(y), [f.coeff for f in factors])
                pre_op = TensorProdSPAMVec("prep", [f.pre_ops[0] for f in factors
                                                    if (f.pre_ops[0] is not None)])
                post_op = TensorProdSPAMVec("prep", [f.post_ops[0] for f in factors
                                                     if (f.post_ops[0] is not None)])
                term = _term.RankOnePolyPrepTerm.simple_init(coeff, pre_op, post_op, self._evotype)

                if not collapsible:  # then may need to add more ops.  Assume factor ops are clifford gates
                    # Embed each factors ops according to their target qubit(s) and just daisy chain them
                    stateSpaceLabels = tuple(range(total_nQ)); curQ = 0
                    for f, nq in zip(factors, fnq):
                        targetLabels = tuple(range(curQ, curQ + nq)); curQ += nq
                        term.pre_ops.extend([_EmbeddedGateMap(stateSpaceLabels, targetLabels, op)
                                             for op in f.pre_ops[1:]])  # embed and add ops
                        term.post_ops.extend([_EmbeddedGateMap(stateSpaceLabels, targetLabels, op)
                                              for op in f.post_ops[1:]])  # embed and add ops

                terms.append(term)

        if return_coeff_polys:
            def _decompose_indices(x):
                return tuple(_modelmember._decompose_gpindices(
                    self.gpindices, _np.array(x, _np.int64)))

            poly_coeffs = [t.coeff.map_indices(_decompose_indices) for t in terms]  # with *local* indices
            tapes = [poly.compact(complex_coeff_tape=True) for poly in poly_coeffs]
            if len(tapes) > 0:
                vtape = _np.concatenate([t[0] for t in tapes])
                ctape = _np.concatenate([t[1] for t in tapes])
            else:
                vtape = _np.empty(0, _np.int64)
                ctape = _np.empty(0, complex)
            coeffs_as_compact_polys = (vtape, ctape)
            #self.local_term_poly_coeffs[order] = coeffs_as_compact_polys #FUTURE?
            return terms, coeffs_as_compact_polys
        else:
            return terms  # Cache terms in FUTURE?

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return self.Np

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        if self._prep_or_effect == "prep":
            return _np.concatenate([fct.to_vector() for fct in self.factors], axis=0)
        else:
            raise ValueError(("'`to_vector` should not be called on effect-like"
                              " TensorProdSPAMVecs (instead it should be called"
                              " on the POVM)"))

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        if self._prep_or_effect == "prep":
            for sv in self.factors:
                sv.from_vector(v[sv.gpindices], close, nodirty)  # factors hold local indices

        elif all([self.effectLbls[i] == list(povm.keys())[0]
                  for i, povm in enumerate(self.factors)]):
            #then this is the *first* vector in the larger TensorProdPOVM
            # and we should initialize all of the factorPOVMs
            for povm in self.factors:
                local_inds = _modelmember._decompose_gpindices(
                    self.gpindices, povm.gpindices)
                povm.from_vector(v[local_inds], close, nodirty)

        #Update representation, which may be a dense matrix or
        # just fast-kron arrays or a stabilizer state.
        self._update_rep()

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.
        An empty 2D array in the StaticSPAMVec case (num_params == 0).

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        assert(self._evotype in ("statevec", "densitymx"))
        typ = complex if self._evotype == "statevec" else 'd'
        derivMx = _np.zeros((self.dim, self.num_params()), typ)

        #Product rule to compute jacobian
        for i, fct in enumerate(self.factors):  # loop over the spamvec/povm we differentiate wrt
            vec = fct if (self._prep_or_effect == "prep") else fct[self.effectLbls[i]]

            if vec.num_params() == 0: continue  # no contribution
            deriv = vec.deriv_wrt_params(None)  # TODO: use filter?? / make relative to this gate...
            deriv.shape = (vec.dim, vec.num_params())

            if i > 0:  # factors before ith
                if self._prep_or_effect == "prep":
                    pre = self.factors[0].todense()
                    for vecA in self.factors[1:i]:
                        pre = _np.kron(pre, vecA.todense())
                else:
                    pre = self.factors[0][self.effectLbls[0]].todense()
                    for j, fctA in enumerate(self.factors[1:i], start=1):
                        pre = _np.kron(pre, fctA[self.effectLbls[j]].todense())
                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
                if self._prep_or_effect == "prep":
                    post = self.factors[i + 1].todense()
                    for vecA in self.factors[i + 2:]:
                        post = _np.kron(post, vecA.todense())
                else:
                    post = self.factors[i + 1][self.effectLbls[i + 1]].todense()
                    for j, fctA in enumerate(self.factors[i + 2:], start=i + 2):
                        post = _np.kron(post, fctA[self.effectLbls[j]].todense())
                deriv = _np.kron(deriv, post[:, None])  # add a dummy 1-dim to 'post' and do kron properly...

            if self._prep_or_effect == "prep":
                local_inds = fct.gpindices  # factor vectors hold local indices
            else:  # in effect case, POVM-factors hold global indices (b/c they're meant to be shareable)
                local_inds = _modelmember._decompose_gpindices(
                    self.gpindices, fct.gpindices)

            assert(local_inds is not None), \
                "Error: gpindices has not been initialized for factor %d - cannot compute derivative!" % i
            derivMx[:, local_inds] += deriv

        derivMx.shape = (self.dim, self.num_params())  # necessary?
        if wrtFilter is None:
            return derivMx
        else:
            return _np.take(derivMx, wrtFilter, axis=1)

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        return False

    def __str__(self):
        s = "Tensor product %s vector with length %d\n" % (self._prep_or_effect, self.dim)
        #ar = self.todense()
        #s += _mt.mx_to_string(ar, width=4, prec=2)

        if self._prep_or_effect == "prep":
            # factors are just other SPAMVecs
            s += " x ".join([_mt.mx_to_string(fct.todense(), width=4, prec=2) for fct in self.factors])
        else:
            # factors are POVMs
            s += " x ".join([_mt.mx_to_string(fct[self.effectLbls[i]].todense(), width=4, prec=2)
                             for i, fct in enumerate(self.factors)])
        return s


class PureStateSPAMVec(SPAMVec):
    """
    Encapsulates a SPAM vector that is a pure state but evolves according to
    one of the density matrix evolution types ("denstiymx", "svterm", and
    "cterm").  It is parameterized by a contained pure-state SPAMVec which
    evolves according to a state vector evolution type ("statevec" or
    "stabilizer").
    """

    def __init__(self, pure_state_vec, evotype='densitymx', dm_basis='pp', typ="prep"):
        """
        Initialize a PureStateSPAMVec object.

        Parameters
        ----------
        pure_state_vec : array_like or SPAMVec
            a 1D numpy array or object representing the pure state.  This object
            sets the parameterization and dimension of this SPAM vector (if
            `pure_state_vec`'s dimension is `d`, then this SPAM vector's
            dimension is `d^2`).  Assumed to be a complex vector in the
            standard computational basis.

        evotype : {'densitymx', 'svterm', 'cterm'}
            The evolution type of this SPAMVec.  Note that the evotype of
            `pure_state_vec` must be compatible with this value.  In particular,
            `pure_state_vec` must have an evotype of `"statevec"` (then allowed
            values are `"densitymx"` and `"svterm"`) or `"stabilizer"` (then
            the only allowed value is `"cterm"`).

        dm_basis : {'std', 'gm', 'pp', 'qt'} or Basis object
            The basis for this SPAM vector - that is, for the *density matrix*
            corresponding to `pure_state_vec`.  Allowed values are Matrix-unit
            (std),  Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt)
            (or a custom basis object).
        """
        if not isinstance(pure_state_vec, SPAMVec):
            pure_state_vec = StaticSPAMVec(SPAMVec.convert_to_vector(pure_state_vec), 'statevec')
        self.pure_state_vec = pure_state_vec
        self.basis = dm_basis  # only used for dense conversion

        pure_evo = pure_state_vec._evotype
        if pure_evo == "statevec":
            if evotype not in ("densitymx", "svterm"):
                raise ValueError(("`evotype` arg must be 'densitymx' or 'svterm'"
                                  " when `pure_state_vec` evotype is 'statevec'"))
        elif pure_evo == "stabilizer":
            if evotype not in ("cterm",):
                raise ValueError(("`evotype` arg must be 'densitymx' or 'svterm'"
                                  " when `pure_state_vec` evotype is 'statevec'"))
        else:
            raise ValueError("`pure_state_vec` evotype must be 'statevec' or 'stabilizer' (not '%s')" % pure_evo)

        dim = self.pure_state_vec.dim**2
        rep = dim  # no representation yet... maybe this should be a dense vector
        # (in "densitymx" evotype case -- use todense)? TODO

        #Create representation
        SPAMVec.__init__(self, rep, evotype, typ)

    def todense(self, scratch=None):
        """
        Return this SPAM vector as a (dense) numpy array.  The memory
        in `scratch` maybe used when it is not-None.
        """
        dmVec_std = _gt.state_to_dmvec(self.pure_state_vec.todense())
        return _bt.change_basis(dmVec_std, 'std', self.basis)

    def get_taylor_order_terms(self, order, max_poly_vars=100, return_coeff_polys=False):
        """
        Get the `order`-th order Taylor-expansion terms of this SPAM vector.

        This function either constructs or returns a cached list of the terms at
        the given order.  Each term is "rank-1", meaning that it is a state
        preparation followed by or POVM effect preceded by actions on a
        density matrix `rho` of the form:

        `rho -> A rho B`

        The coefficients of these terms are typically polynomials of the
        SPAMVec's parameters, where the polynomial's variable indices index the
        *global* parameters of the SPAMVec's parent (usually a :class:`Model`)
        , not the SPAMVec's local parameter array (i.e. that returned from
        `to_vector`).


        Parameters
        ----------
        order : int
            The order of terms to get.

        return_coeff_polys : bool
            Whether a parallel list of locally-indexed (using variable indices
            corresponding to *this* object's parameters rather than its parent's)
            polynomial coefficients should be returned as well.

        Returns
        -------
        terms : list
            A list of :class:`RankOneTerm` objects.

        coefficients : list
            Only present when `return_coeff_polys == True`.
            A list of *compact* polynomial objects, meaning that each element
            is a `(vtape,ctape)` 2-tuple formed by concatenating together the
            output of :method:`Polynomial.compact`.
        """
        if self.num_params() > 0:
            raise ValueError(("PureStateSPAMVec.get_taylor_order_terms(...) is only "
                              "implemented for the case when its underlying "
                              "pure state vector has 0 parameters (is static)"))

        if order == 0:  # only 0-th order term exists (assumes static pure_state_vec)
            purevec = self.pure_state_vec
            coeff = _Polynomial({(): 1.0}, max_poly_vars)
            if self._prep_or_effect == "prep":
                terms = [_term.RankOnePolyPrepTerm.simple_init(coeff, purevec, purevec, self._evotype)]
            else:
                terms = [_term.RankOnePolyEffectTerm.simple_init(coeff, purevec, purevec, self._evotype)]

            if return_coeff_polys:
                coeffs_as_compact_polys = coeff.compact(complex_coeff_tape=True)
                return terms, coeffs_as_compact_polys
            else:
                return terms
        else:
            if return_coeff_polys:
                vtape = _np.empty(0, _np.int64)
                ctape = _np.empty(0, complex)
                return [], (vtape, ctape)
            else:
                return []

    #TODO REMOVE
    #def get_direct_order_terms(self, order, base_order):
    #    """
    #    Parameters
    #    ----------
    #    order : int
    #        The order of terms to get.
    #
    #    Returns
    #    -------
    #    list
    #        A list of :class:`RankOneTerm` objects.
    #    """
    #    if self.num_params() > 0:
    #        raise ValueError(("PureStateSPAMVec.get_taylor_order_terms(...) is only "
    #                          "implemented for the case when its underlying "
    #                          "pure state vector has 0 parameters (is static)"))
    #
    #    if order == 0: # only 0-th order term exists (assumes static pure_state_vec)
    #        if self._evotype == "svterm":  tt = "dense"
    #        elif self._evotype == "cterm": tt = "clifford"
    #        else: raise ValueError("Invalid evolution type %s for calling `get_taylor_order_terms`" % self._evotype)
    #
    #        purevec = self.pure_state_vec
    #        terms = [ _term.RankOneTerm(1.0, purevec, purevec, tt) ]
    #    else:
    #        terms = []
    #    return terms

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return self.pure_state_vec.num_params()

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        return self.pure_state_vec.to_vector()

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        self.pure_state_vec.from_vector(v, close, nodirty)
        #Update dense rep if one is created (TODO)

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.
        An empty 2D array in the StaticSPAMVec case (num_params == 0).

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        raise NotImplementedError("Still need to work out derivative calculation of PureStateSPAMVec")

    def has_nonzero_hessian(self):
        """
        Returns whether this SPAM vector has a non-zero Hessian with
        respect to its parameters, i.e. whether it only depends
        linearly on its parameters or not.

        Returns
        -------
        bool
        """
        return self.pure_state_vec.has_nonzero_hessian()

    def __str__(self):
        s = "Pure-state spam vector with length %d holding:\n" % self.dim
        s += "  " + str(self.pure_state_vec)
        return s


class LindbladSPAMVec(SPAMVec):
    """ A Lindblad-parameterized SPAMVec (that is also expandable into terms)"""

    @classmethod
    def from_spamvec_obj(cls, spamvec, typ, paramType="GLND", purevec=None,
                         proj_basis="pp", mxBasis="pp", truncate=True,
                         lazy=False):
        """
        Creates a LindbladSPAMVec from an existing SPAMVec object
        and some additional information.

        This function is different from `from_spam_vector` in that it assumes
        that `spamvec` is a :class:`SPAMVec`-derived object, and if `lazy=True`
        and if `spamvec` is already a matching LindbladSPAMVec, it
        is returned directly.  This routine is primarily used in spam vector
        conversion functions, where conversion is desired only when necessary.

        Parameters
        ----------
        spamvec : SPAMVec
            The spam vector object to "convert" to a
            `LindbladSPAMVec`.

        typ : {"prep","effect"}
            Whether this is a state preparation or POVM effect vector.

        paramType : str, optional
            The high-level "parameter type" of the gate to create.  This
            specifies both which Lindblad parameters are included and what
            type of evolution is used.  Examples of valid values are
            `"CPTP"`, `"H+S"`, `"S terms"`, and `"GLND clifford terms"`.

        purevec : numpy array or SPAMVec object, optional
            A SPAM vector which represents a pure-state, taken as the "ideal"
            reference state when constructing the error generator of the
            returned `LindbladSPAMVec`.  Note that this vector
            still acts on density matrices (if it's a SPAMVec it should have
            a "densitymx", "svterm", or "cterm" evolution type, and if it's
            a numpy array it should have the same dimension as `spamvec`).
            If None, then it is taken to be `spamvec`, and so `spamvec` must
            represent a pure state in this case.

        proj_basis: {'std', 'gm', 'pp', 'qt'}, list of matrices, or Basis object
            The basis used to construct the Lindblad-term error generators onto
            which the SPAM vector's error generator is projected.  Allowed values
            are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
            and Qutrit (qt), list of numpy arrays, or a custom basis object.

        mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
            The source and destination basis, respectively.  Allowed
            values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
            and Qutrit (qt) (or a custom basis object).

        truncate : bool, optional
            Whether to truncate the projections onto the Lindblad terms in
            order to meet constraints (e.g. to preserve CPTP) when necessary.
            If False, then an error is thrown when the given `spamvec` cannot
            be realized by the specified set of Lindblad projections.

        lazy : bool, optional
            If True, then if `spamvec` is already a LindbladSPAMVec
            with the requested details (given by the other arguments), then
            `spamvec` is returned directly and no conversion/copying is
            performed. If False, then a new object is always returned.

        Returns
        -------
        LindbladSPAMVec
        """

        if not isinstance(spamvec, SPAMVec):
            spamvec = StaticSPAMVec(spamvec, typ=typ)  # assume spamvec is just a vector

        if purevec is None:
            purevec = spamvec  # right now, we don't try to extract a "closest pure vec"
            # to spamvec - below will fail if spamvec isn't pure.
        elif not isinstance(purevec, SPAMVec):
            purevec = StaticSPAMVec(purevec, typ=typ)  # assume spamvec is just a vector

        #Break paramType in to a "base" type and an evotype
        from .operation import LindbladOp as _LPGMap
        bTyp, evotype, nonham_mode, param_mode = _LPGMap.decomp_paramtype(paramType)

        ham_basis = proj_basis if (("H+" in bTyp) or bTyp in ("CPTP", "GLND")) else None
        nonham_basis = proj_basis

        def beq(b1, b2):
            """ Check if bases have equal names """
            b1 = b1.name if isinstance(b1, _Basis) else b1
            b2 = b2.name if isinstance(b2, _Basis) else b2
            return b1 == b2

        def normeq(a, b):
            if a is None and b is None: return True
            if a is None or b is None: return False
            return _mt.safenorm(a - b) < 1e-6  # what about possibility of Clifford gates?

        if isinstance(spamvec, LindbladSPAMVec) \
           and spamvec._evotype == evotype and spamvec.typ == typ \
           and beq(ham_basis, spamvec.error_map.ham_basis) and beq(nonham_basis, spamvec.error_map.other_basis) \
           and param_mode == spamvec.error_map.param_mode and nonham_mode == spamvec.error_map.nonham_mode \
           and beq(mxBasis, spamvec.error_map.matrix_basis) and lazy:
            #normeq(gate.pure_state_vec,purevec) \ # TODO: more checks for equality?!
            return spamvec  # no creation necessary!
        else:
            #Convert vectors (if possible) to SPAMVecs
            # of the appropriate evotype and 0 params.
            bDiff = spamvec is not purevec
            spamvec = _convert_to_lindblad_base(spamvec, typ, evotype, mxBasis)
            purevec = _convert_to_lindblad_base(purevec, typ, evotype, mxBasis) if bDiff else spamvec
            assert(spamvec._evotype == evotype)
            assert(purevec._evotype == evotype)

            return cls.from_spam_vector(
                spamvec, purevec, typ, ham_basis, nonham_basis,
                param_mode, nonham_mode, truncate, mxBasis, evotype)

    @classmethod
    def from_spam_vector(cls, spamVec, pureVec, typ,
                         ham_basis="pp", nonham_basis="pp", param_mode="cptp",
                         nonham_mode="all", truncate=True, mxBasis="pp",
                         evotype="densitymx"):
        """
        Creates a Lindblad-parameterized spamvec from a state vector and a basis
        which specifies how to decompose (project) the vector's error generator.

        spamVec : SPAMVec
            the SPAM vector to initialize from.  The error generator that
            tranforms `pureVec` into `spamVec` forms the parameterization
            of the returned LindbladSPAMVec.

        pureVec : numpy array or SPAMVec
            An array or SPAMVec in the *full* density-matrix space (this
            vector will have the same dimension as `spamVec` - 4 in the case
            of a single qubit) which represents a pure-state preparation or
            projection.  This is used as the "base" preparation/projection
            when computing the error generator that will be parameterized.
            Note that this argument must be specified, as there is no natural
            default value (like the identity in the case of gates).

        typ : {"prep","effect"}
            Whether this is a state preparation or POVM effect vector.

        ham_basis: {'std', 'gm', 'pp', 'qt'}, list of matrices, or Basis object
            The basis is used to construct the Hamiltonian-type lindblad error
            Allowed values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
            and Qutrit (qt), list of numpy arrays, or a custom basis object.

        nonham_basis: {'std', 'gm', 'pp', 'qt'}, list of matrices, or Basis object
            The basis is used to construct the Stochastic-type lindblad error
            Allowed values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
            and Qutrit (qt), list of numpy arrays, or a custom basis object.

        param_mode : {"unconstrained", "cptp", "depol", "reldepol"}
            Describes how the Lindblad coefficients/projections relate to the
            SPAM vector's parameter values.  Allowed values are:
            `"unconstrained"` (coeffs are independent unconstrained parameters),
            `"cptp"` (independent parameters but constrained so map is CPTP),
            `"reldepol"` (all non-Ham. diagonal coeffs take the *same* value),
            `"depol"` (same as `"reldepol"` but coeffs must be *positive*)

        nonham_mode : {"diagonal", "diag_affine", "all"}
            Which non-Hamiltonian Lindblad projections are potentially non-zero.
            Allowed values are: `"diagonal"` (only the diagonal Lind. coeffs.),
            `"diag_affine"` (diagonal coefficients + affine projections), and
            `"all"` (the entire matrix of coefficients is allowed).

        truncate : bool, optional
            Whether to truncate the projections onto the Lindblad terms in
            order to meet constraints (e.g. to preserve CPTP) when necessary.
            If False, then an error is thrown when the given `gate` cannot
            be realized by the specified set of Lindblad projections.

        mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
            The source and destination basis, respectively.  Allowed
            values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
            and Qutrit (qt) (or a custom basis object).

        evotype : {"densitymx","svterm","cterm"}
            The evolution type of the spamvec being constructed.  `"densitymx"` is
            usual Lioville density-matrix-vector propagation via matrix-vector
            products.  `"svterm"` denotes state-vector term-based evolution
            (spamvec is obtained by evaluating the rank-1 terms up to
            some order).  `"cterm"` is similar but stabilizer states.

        Returns
        -------
        LindbladSPAMVec
        """
        #Compute a (errgen, pureVec) pair from the given
        # (spamVec, pureVec) pair.

        assert(pureVec is not None), "Must supply `pureVec`!"  # since there's no good default?

        if not isinstance(spamVec, SPAMVec):
            spamVec = StaticSPAMVec(spamVec, evotype, typ)  # assume spamvec is just a vector
        if not isinstance(pureVec, SPAMVec):
            pureVec = StaticSPAMVec(pureVec, evotype, typ)  # assume spamvec is just a vector
        d2 = pureVec.dim

        #Determine whether we're using sparse bases or not
        sparse = None
        if ham_basis is not None:
            if isinstance(ham_basis, _Basis): sparse = ham_basis.sparse
            elif not isinstance(ham_basis, str) and len(ham_basis) > 0:
                sparse = _sps.issparse(ham_basis[0])
        if sparse is None and nonham_basis is not None:
            if isinstance(nonham_basis, _Basis): sparse = nonham_basis.sparse
            elif not isinstance(nonham_basis, str) and len(nonham_basis) > 0:
                sparse = _sps.issparse(nonham_basis[0])
        if sparse is None: sparse = False  # the default

        if spamVec is None or spamVec is pureVec:
            if sparse: errgen = _sps.csr_matrix((d2, d2), dtype='d')
            else: errgen = _np.zeros((d2, d2), 'd')
        else:
            #Construct "spam error generator" by comparing *dense* vectors
            pvdense = pureVec.todense()
            svdense = spamVec.todense()
            errgen = _gt.spam_error_generator(svdense, pvdense, mxBasis)
            if sparse: errgen = _sps.csr_matrix(errgen)

        assert(pureVec._evotype == evotype), "`pureVec` must have evotype == '%s'" % evotype

        from .operation import LindbladErrorgen as _LErrorgen
        from .operation import LindbladOp as _LPGMap
        from .operation import LindbladDenseOp as _LPOp

        errgen = _LErrorgen.from_error_generator(errgen, ham_basis,
                                                 nonham_basis, param_mode, nonham_mode,
                                                 mxBasis, truncate, evotype)
        errcls = _LPOp if (pureVec.dim <= 64 and evotype == "densitymx") else _LPGMap
        errmap = errcls(None, errgen)

        return cls(pureVec, errmap, typ)

    @classmethod
    def from_lindblad_terms(cls, pureVec, Ltermdict, typ, basisdict=None,
                            param_mode="cptp", nonham_mode="all", truncate=True,
                            mxBasis="pp", evotype="densitymx"):
        """
        Create a Lindblad-parameterized spamvec with a given set of Lindblad terms.

        Parameters
        ----------
        pureVec : numpy array or SPAMVec
            An array or SPAMVec in the *full* density-matrix space (this
            vector will have dimension 4 in the case of a single qubit) which
            represents a pure-state preparation or projection.  This is used as
            the "base" preparation or projection that is followed or preceded
            by, respectively, the parameterized Lindblad-form error generator.

        Ltermdict : dict
            A dictionary specifying which Linblad terms are present in the gate
            parameteriztion.  Keys are `(termType, basisLabel1, <basisLabel2>)`
            tuples, where `termType` can be `"H"` (Hamiltonian), `"S"`
            (Stochastic), or `"A"` (Affine).  Hamiltonian and Affine terms always
            have a single basis label (so key is a 2-tuple) whereas Stochastic
            tuples with 1 basis label indicate a *diagonal* term, and are the
            only types of terms allowed when `nonham_mode != "all"`.  Otherwise,
            Stochastic term tuples can include 2 basis labels to specify
            "off-diagonal" non-Hamiltonian Lindblad terms.  Basis labels can be
            strings or integers.  Values are complex coefficients (error rates).

        typ : {"prep","effect"}
            Whether this is a state preparation or POVM effect vector.

        basisdict : dict, optional
            A dictionary mapping the basis labels (strings or ints) used in the
            keys of `Ltermdict` to basis matrices (numpy arrays or Scipy sparse
            matrices).

        param_mode : {"unconstrained", "cptp", "depol", "reldepol"}
            Describes how the Lindblad coefficients/projections relate to the
            SPAM vector's parameter values.  Allowed values are:
            `"unconstrained"` (coeffs are independent unconstrained parameters),
            `"cptp"` (independent parameters but constrained so map is CPTP),
            `"reldepol"` (all non-Ham. diagonal coeffs take the *same* value),
            `"depol"` (same as `"reldepol"` but coeffs must be *positive*)

        nonham_mode : {"diagonal", "diag_affine", "all"}
            Which non-Hamiltonian Lindblad projections are potentially non-zero.
            Allowed values are: `"diagonal"` (only the diagonal Lind. coeffs.),
            `"diag_affine"` (diagonal coefficients + affine projections), and
            `"all"` (the entire matrix of coefficients is allowed).

        truncate : bool, optional
            Whether to truncate the projections onto the Lindblad terms in
            order to meet constraints (e.g. to preserve CPTP) when necessary.
            If False, then an error is thrown when the given dictionary of
            Lindblad terms doesn't conform to the constrains.

        mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
            The source and destination basis, respectively.  Allowed
            values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
            and Qutrit (qt) (or a custom basis object).

        evotype : {"densitymx","svterm","cterm"}
            The evolution type of the spamvec being constructed.  `"densitymx"` is
            usual Lioville density-matrix-vector propagation via matrix-vector
            products.  `"svterm"` denotes state-vector term-based evolution
            (spamvec is obtained by evaluating the rank-1 terms up to
            some order).  `"cterm"` is similar but stabilizer states.

        Returns
        -------
        LindbladOp
        """
        #Need a dimension for error map construction (basisdict could be completely empty)
        if not isinstance(pureVec, SPAMVec):
            pureVec = StaticSPAMVec(pureVec, evotype, typ)  # assume spamvec is just a vector
        d2 = pureVec.dim

        from .operation import LindbladOp as _LPGMap
        errmap = _LPGMap(d2, Ltermdict, basisdict, param_mode, nonham_mode,
                         truncate, mxBasis, evotype)
        return cls(pureVec, errmap, typ)

    def __init__(self, pureVec, errormap, typ):
        """
        Initialize a LindbladSPAMVec object.

        Essentially a pure state preparation or projection that is followed
        or preceded by, respectively, the action of LindbladDenseOp.

        Parameters
        ----------
        pureVec : numpy array or SPAMVec
            An array or SPAMVec in the *full* density-matrix space (this
            vector will have dimension 4 in the case of a single qubit) which
            represents a pure-state preparation or projection.  This is used as
            the "base" preparation or projection that is followed or preceded
            by, respectively, the parameterized Lindblad-form error generator.
            (This argument is *not* copied if it is a SPAMVec.  A numpy array
             is converted to a new StaticSPAMVec.)

        errormap : MapOperator
            The error generator action and parameterization, encapsulated in
            a gate object.  Usually a :class:`LindbladOp`
            or :class:`ComposedOp` object.  (This argument is *not* copied,
            to allow LindbladSPAMVecs to share error generator
            parameters with other gates and spam vectors.)

        typ : {"prep","effect"}
            Whether this is a state preparation or POVM effect vector.
        """
        from .operation import LindbladOp as _LPGMap
        evotype = errormap._evotype
        assert(evotype in ("densitymx", "svterm", "cterm")), \
            "Invalid evotype: %s for %s" % (evotype, self.__class__.__name__)

        if not isinstance(pureVec, SPAMVec):
            pureVec = StaticSPAMVec(pureVec, evotype, typ)  # assume spamvec is just a vector

        assert(pureVec._evotype == evotype), \
            "`pureVec` evotype must match `errormap` ('%s' != '%s')" % (pureVec._evotype, evotype)
        assert(pureVec.num_params() == 0), "`pureVec` 'reference' must have *zero* parameters!"

        d2 = pureVec.dim
        self.state_vec = pureVec
        self.error_map = errormap
        self.terms = {} if evotype in ("svterm", "cterm") else None
        self.local_term_poly_coeffs = {} if evotype in ("svterm", "cterm") else None
        # TODO REMOVE self.direct_terms = {} if evotype in ("svterm","cterm") else None
        # TODO REMOVE self.direct_term_poly_coeffs = {} if evotype in ("svterm","cterm") else None

        #Create representation
        if evotype == "densitymx":
            assert(self.state_vec._prep_or_effect == typ), "LindbladSPAMVec prep/effect mismatch with given statevec!"

            if typ == "prep":
                dmRep = self.state_vec._rep
                errmapRep = self.error_map._rep
                rep = errmapRep.acton(dmRep)  # FUTURE: do this acton in place somehow? (like C-reps do)
                #maybe make a special _Errgen *state* rep?

            else:  # effect
                dmEffectRep = self.state_vec._rep
                errmapRep = self.error_map._rep
                rep = replib.DMEffectRep_Errgen(errmapRep, dmEffectRep, id(self.error_map))
                # an effect that applies a *named* errormap before computing with dmEffectRep

        else:
            rep = d2  # no representations for term-based evotypes

        SPAMVec.__init__(self, rep, evotype, typ)  # sets self.dim

    def _update_rep(self):
        if self._evotype == "densitymx":
            if self._prep_or_effect == "prep":
                # _rep is a DMStateRep
                dmRep = self.state_vec._rep
                errmapRep = self.error_map._rep
                self._rep.base[:] = errmapRep.acton(dmRep).base[:]  # copy from "new_rep"

            else:  # effect
                # _rep is a DMEffectRep_Errgen, which just holds references to the
                # effect and error map's representations (which we assume have been updated)
                # so there's no need to update anything here
                pass

    def submembers(self):
        """
        Get the ModelMember-derived objects contained in this one.

        Returns
        -------
        list
        """
        return [self.error_map]

    def copy(self, parent=None):
        """
        Copy this object.

        Returns
        -------
        LinearOperator
            A copy of this object.
        """
        # We need to override this method so that embedded gate has its
        # parent reset correctly.
        cls = self.__class__  # so that this method works for derived classes too
        copyOfMe = cls(self.state_vec, self.error_map.copy(parent), self._prep_or_effect)
        return self._copy_gpindices(copyOfMe, parent)

    def set_gpindices(self, gpindices, parent, memo=None):
        """
        Set the parent and indices into the parent's parameter vector that
        are used by this ModelMember object.

        Parameters
        ----------
        gpindices : slice or integer ndarray
            The indices of this objects parameters in its parent's array.

        parent : Model or ModelMember
            The parent whose parameter array gpindices references.

        Returns
        -------
        None
        """
        self.terms = {}  # clear terms cache since param indices have changed now
        self.local_term_poly_coeffs = {}
        # TODO REMOVE self.direct_terms = {}
        # TODO REMOVE self.direct_term_poly_coeffs = {}
        _modelmember.ModelMember.set_gpindices(self, gpindices, parent, memo)

    def todense(self, scratch=None):
        """
        Return this SPAM vector as a (dense) numpy array.  The memory
        in `scratch` maybe used when it is not-None.
        """
        if self._prep_or_effect == "prep":
            #error map acts on dmVec
            return _np.dot(self.error_map.todense(), self.state_vec.todense())
        else:
            #Note: if this is an effect vector, self.error_map is the
            # map that acts on the *state* vector before dmVec acts
            # as an effect:  E.T -> dot(E.T,errmap) ==> E -> dot(errmap.T,E)
            return _np.dot(self.error_map.todense().conjugate().T, self.state_vec.todense())

    #def torep(self, typ, outvec=None):
    #    """
    #    Return a "representation" object for this SPAMVec.
    #
    #    Such objects are primarily used internally by pyGSTi to compute
    #    things like probabilities more efficiently.
    #
    #    Returns
    #    -------
    #    StateRep
    #    """
    #    if self._evotype == "densitymx":
    #
    #        if typ == "prep":
    #            dmRep = self.state_vec.torep(typ)
    #            errmapRep = self.error_map.torep()
    #            return errmapRep.acton(dmRep)  # FUTURE: do this acton in place somehow? (like C-reps do)
    #            #maybe make a special _Errgen *state* rep?
    #
    #        else:  # effect
    #            dmEffectRep = self.state_vec.torep(typ)
    #            errmapRep = self.error_map.torep()
    #            return replib.DMEffectRep_Errgen(errmapRep, dmEffectRep, id(self.error_map))
    #            # an effect that applies a *named* errormap before computing with dmEffectRep
    #
    #    else:
    #        #framework should not be calling "torep" on states w/a term-based evotype...
    #        # they should call torep on the *terms* given by get_taylor_order_terms(...)
    #        raise ValueError("Invalid evotype '%s' for %s.torep(...)" %
    #                         (self._evotype, self.__class__.__name__))

    def get_taylor_order_terms(self, order, max_poly_vars=100, return_coeff_polys=False):
        """
        Get the `order`-th order Taylor-expansion terms of this SPAM vector.

        This function either constructs or returns a cached list of the terms at
        the given order.  Each term is "rank-1", meaning that it is a state
        preparation followed by or POVM effect preceded by actions on a
        density matrix `rho` of the form:

        `rho -> A rho B`

        The coefficients of these terms are typically polynomials of the
        SPAMVec's parameters, where the polynomial's variable indices index the
        *global* parameters of the SPAMVec's parent (usually a :class:`Model`)
        , not the SPAMVec's local parameter array (i.e. that returned from
        `to_vector`).


        Parameters
        ----------
        order : int
            The order of terms to get.

        return_coeff_polys : bool
            Whether a parallel list of locally-indexed (using variable indices
            corresponding to *this* object's parameters rather than its parent's)
            polynomial coefficients should be returned as well.

        Returns
        -------
        terms : list
            A list of :class:`RankOneTerm` objects.

        coefficients : list
            Only present when `return_coeff_polys == True`.
            A list of *compact* polynomial objects, meaning that each element
            is a `(vtape,ctape)` 2-tuple formed by concatenating together the
            output of :method:`Polynomial.compact`.
        """
        if order not in self.terms:
            if self._evotype not in ('svterm', 'cterm'):
                raise ValueError("Invalid evolution type %s for calling `get_taylor_order_terms`" % self._evotype)
            assert(self.gpindices is not None), "LindbladSPAMVec must be added to a Model before use!"

            state_terms = self.state_vec.get_taylor_order_terms(0, max_poly_vars); assert(len(state_terms) == 1)
            stateTerm = state_terms[0]
            err_terms, cpolys = self.error_map.get_taylor_order_terms(order, max_poly_vars, True)
            if self._prep_or_effect == "prep":
                terms = [_term.compose_terms((stateTerm, t)) for t in err_terms]  # t ops occur *after* stateTerm's
            else:  # "effect"
                # Effect terms are special in that all their pre/post ops act in order on the *state* before the final
                # effect is used to compute a probability.  Thus, constructing the same "terms" as above works here
                # too - the difference comes when this SPAMVec is used as an effect rather than a prep.
                terms = [_term.compose_terms((stateTerm, t)) for t in err_terms]  # t ops occur *after* stateTerm's

            #OLD: now this is done within calculator when possible b/c not all terms can be collapsed
            #terms = [ t.collapse() for t in terms ] # collapse terms for speed
            # - resulting in terms with just a single pre/post op, each == a pure state

            #assert(stateTerm.coeff == Polynomial_1.0) # TODO... so can assume local polys are same as for errorgen
            self.local_term_poly_coeffs[order] = cpolys
            self.terms[order] = terms

        if return_coeff_polys:
            return self.terms[order], self.local_term_poly_coeffs[order]
        else:
            return self.terms[order]

    def get_taylor_order_terms_above_mag(self, order, max_poly_vars, min_term_mag):
        state_terms = self.state_vec.get_taylor_order_terms(0, max_poly_vars); assert(len(state_terms) == 1)
        stateTerm = state_terms[0]
        stateTerm = stateTerm.copy_with_magnitude(1.0)
        #assert(stateTerm.coeff == Polynomial_1.0) # TODO... so can assume local polys are same as for errorgen

        err_terms = self.error_map.get_taylor_order_terms_above_mag(
            order, max_poly_vars, min_term_mag / stateTerm.magnitude)

        #This gives the appropriate logic, but *both* prep or effect results in *same* expression, so just run it:
        #if self._prep_or_effect == "prep":
        #    terms = [_term.compose_terms((stateTerm, t)) for t in err_terms]  # t ops occur *after* stateTerm's
        #else:  # "effect"
        #    # Effect terms are special in that all their pre/post ops act in order on the *state* before the final
        #    # effect is used to compute a probability.  Thus, constructing the same "terms" as above works here
        #    # too - the difference comes when this SPAMVec is used as an effect rather than a prep.
        #    terms = [_term.compose_terms((stateTerm, t)) for t in err_terms]  # t ops occur *after* stateTerm's
        terms = [_term.compose_terms_with_mag((stateTerm, t), stateTerm.magnitude * t.magnitude)
                 for t in err_terms]  # t ops occur *after* stateTerm's
        return terms

    def get_total_term_magnitude(self):
        """
        Get the total (sum) of the magnitudes of all this SPAM vector's terms.

        The magnitude of a term is the absolute value of its coefficient, so
        this function returns the number you'd get from summing up the
        absolute-coefficients of all the Taylor terms (at all orders!) you
        get from expanding this SPAM vector in a Taylor series.

        Returns
        -------
        float
        """
        # return (sum of absvals of *all* term coeffs)
        return self.error_map.get_total_term_magnitude()  # error map is only part with terms

    def get_total_term_magnitude_deriv(self):
        """
        Get the derivative of the total (sum) of the magnitudes of all this
        operator's terms with respect to the operators (local) parameters.

        Returns
        -------
        numpy array
            An array of length self.num_params()
        """
        return self.error_map.get_total_term_magnitude_deriv()

    def deriv_wrt_params(self, wrtFilter=None):
        """
        Construct a matrix whose columns are the derivatives of the SPAM vector
        with respect to a single param.  Thus, each column is of length
        get_dimension and there is one column per SPAM vector parameter.

        Returns
        -------
        numpy array
            Array of derivatives, shape == (dimension, num_params)
        """
        dmVec = self.state_vec.todense()

        derrgen = self.error_map.deriv_wrt_params(wrtFilter)  # shape (dim*dim, nParams)
        derrgen.shape = (self.dim, self.dim, derrgen.shape[1])  # => (dim,dim,nParams)

        if self._prep_or_effect == "prep":
            #derror map acts on dmVec
            #return _np.einsum("ijk,j->ik", derrgen, dmVec) # return shape = (dim,nParams)
            return _np.tensordot(derrgen, dmVec, (1, 0))  # return shape = (dim,nParams)
        else:
            # self.error_map acts on the *state* vector before dmVec acts
            # as an effect:  E.dag -> dot(E.dag,errmap) ==> E -> dot(errmap.dag,E)
            #return _np.einsum("jik,j->ik", derrgen.conjugate(), dmVec) # return shape = (dim,nParams)
            return _np.tensordot(derrgen.conjugate(), dmVec, (0, 0))  # return shape = (dim,nParams)

    def hessian_wrt_params(self, wrtFilter1=None, wrtFilter2=None):
        """
        Construct the Hessian of this SPAM vector with respect to its parameters.

        This function returns a tensor whose first axis corresponds to the
        flattened operation matrix and whose 2nd and 3rd axes correspond to the
        parameters that are differentiated with respect to.

        Parameters
        ----------
        wrtFilter1, wrtFilter2 : list
            Lists of indices of the paramters to take first and second
            derivatives with respect to.  If None, then derivatives are
            taken with respect to all of the vectors's parameters.

        Returns
        -------
        numpy array
            Hessian with shape (dimension, num_params1, num_params2)
        """
        dmVec = self.state_vec.todense()

        herrgen = self.error_map.hessian_wrt_params(wrtFilter1, wrtFilter2)  # shape (dim*dim, nParams1, nParams2)
        herrgen.shape = (self.dim, self.dim, herrgen.shape[1], herrgen.shape[2])  # => (dim,dim,nParams1, nParams2)

        if self._prep_or_effect == "prep":
            #derror map acts on dmVec
            #return _np.einsum("ijkl,j->ikl", herrgen, dmVec) # return shape = (dim,nParams)
            return _np.tensordot(herrgen, dmVec, (1, 0))  # return shape = (dim,nParams)
        else:
            # self.error_map acts on the *state* vector before dmVec acts
            # as an effect:  E.dag -> dot(E.dag,errmap) ==> E -> dot(errmap.dag,E)
            #return _np.einsum("jikl,j->ikl", herrgen.conjugate(), dmVec) # return shape = (dim,nParams)
            return _np.tensordot(herrgen.conjugate(), dmVec, (0, 0))  # return shape = (dim,nParams)

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return self.error_map.num_params()

    def to_vector(self):
        """
        Extract a vector of the underlying gate parameters from this gate.

        Returns
        -------
        numpy array
            a 1D numpy array with length == num_params().
        """
        return self.error_map.to_vector()

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the gate using a vector of its parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params().

        Returns
        -------
        None
        """
        self.error_map.from_vector(v, close, nodirty)
        self._update_rep()
        if not nodirty: self.dirty = True

    def transform(self, S, typ):
        """
        Update SPAM (column) vector V as inv(S) * V or S^T * V for preparation
        or  effect SPAM vectors, respectively.

        Note that this is equivalent to state preparation vectors getting
        mapped: `rho -> inv(S) * rho` and the *transpose* of effect vectors
        being mapped as `E^T -> E^T * S`.

        Generally, the transform function updates the *parameters* of
        the SPAM vector such that the resulting vector is altered as
        described above.  If such an update cannot be done (because
        the gate parameters do not allow for it), ValueError is raised.

        Parameters
        ----------
        S : GaugeGroupElement
            A gauge group element which specifies the "S" matrix
            (and it's inverse) used in the above similarity transform.

        typ : { 'prep', 'effect' }
            Which type of SPAM vector is being transformed (see above).
        """
        #Defer everything to LindbladOp's
        # `spam_tranform` function, which applies either
        # error_map -> inv(S) * error_map ("prep" case) OR
        # error_map -> error_map * S      ("effect" case)
        self.error_map.spam_transform(S, typ)
        self._update_rep()
        self.dirty = True

    def depolarize(self, amount):
        """
        Depolarize this SPAM vector by the given `amount`.

        Generally, the depolarize function updates the *parameters* of
        the SPAMVec such that the resulting vector is depolarized.  If
        such an update cannot be done (because the gate parameters do not
        allow for it), ValueError is raised.

        Parameters
        ----------
        amount : float or tuple
            The amount to depolarize by.  If a tuple, it must have length
            equal to one less than the dimension of the spam vector. All but
            the first element of the spam vector (often corresponding to the
            identity element) are multiplied by `amount` (if a float) or
            the corresponding `amount[i]` (if a tuple).

        Returns
        -------
        None
        """
        self.error_map.depolarize(amount)
        self._update_rep()


class StabilizerSPAMVec(SPAMVec):
    """
    A stabilizer state preparation represented internally using a compact
    representation of its stabilizer group.
    """

    @classmethod
    def from_dense_purevec(cls, purevec):
        """
        Create a new StabilizerSPAMVec from a pure-state vector.

        Currently, purevec must be a single computational basis state (it
        cannot be a superpostion of multiple of them).

        Parameters
        ----------
        purevec : numpy.ndarray
            A complex-valued state vector specifying a pure state in the
            standard computational basis.  This vector has length 2^n for
            n qubits.

        Returns
        -------
        StabilizerSPAMVec
        """
        nqubits = int(round(_np.log2(len(purevec))))
        v = (_np.array([1, 0], 'd'), _np.array([0, 1], 'd'))  # (v0,v1)
        for zvals in _itertools.product(*([(0, 1)] * nqubits)):
            testvec = _functools.reduce(_np.kron, [v[i] for i in zvals])
            if _np.allclose(testvec, purevec.flat):
                return cls(nqubits, zvals)
        raise ValueError(("Given `purevec` must be a z-basis product state - "
                          "cannot construct StabilizerSPAMVec"))

    def __init__(self, nqubits, zvals=None, sframe=None):
        """
        Initialize a StabilizerSPAMVec object.

        Parameters
        ----------
        nqubits : int
            Number of qubits

        zvals : iterable, optional
            An iterable over anything that can be cast as True/False
            to indicate the 0/1 value of each qubit in the Z basis.
            If None, the all-zeros state is created.

        sframe : StabilizerFrame, optional
            A complete stabilizer frame to initialize this state from.
            If this is not None, then `nqubits` and `zvals` must be None.
        """
        if sframe is not None:
            assert(nqubits is None and zvals is None), "`nqubits` and `zvals` must be None when `sframe` isn't!"
            self.sframe = sframe
        else:
            self.sframe = _stabilizer.StabilizerFrame.from_zvals(nqubits, zvals)
        rep = self.sframe.torep()  # dim == 2**nqubits
        SPAMVec.__init__(self, rep, "stabilizer", "prep")

    def todense(self, scratch=None):
        """
        Return this SPAM vector as a (dense) numpy array of shape
        (2^(nqubits), 1).  The memory in `scratch` maybe used when
        it is not-None.
        """
        statevec = self.sframe.to_statevec()
        statevec.shape = (statevec.size, 1)
        return statevec

    #def torep(self, typ, outvec=None):
    #    """
    #    Return a "representation" object for this SPAMVec.
    #
    #    Such objects are primarily used internally by pyGSTi to compute
    #    things like probabilities more efficiently.
    #
    #    Returns
    #    -------
    #    SBStateRep
    #    """
    #    return self.sframe.torep()

    def __str__(self):
        s = "Stabilizer spam vector for %d qubits with rep:\n" % (self.sframe.nqubits)
        s += str(self.sframe)
        return s


class StabilizerEffectVec(SPAMVec):  # FUTURE: merge this with ComptationalSPAMVec (w/evotype == "stabilizer")?
    """
    A dummy SPAM vector that points to a set/product of 1-qubit POVM
    outcomes from stabilizer-state measurements.
    """

    @classmethod
    def from_dense_purevec(cls, purevec):
        """
        Create a new StabilizerEffectVec from a pure-state vector.

        Currently, purevec must be a single computational basis state (it
        cannot be a superpostion of multiple of them).

        Parameters
        ----------
        purevec : numpy.ndarray
            A complex-valued state vector specifying a pure state in the
            standard computational basis.  This vector has length 2^n for
            n qubits.

        Returns
        -------
        StabilizerSPAMVec
        """
        nqubits = int(round(_np.log2(len(purevec))))
        v = (_np.array([1, 0], 'd'), _np.array([0, 1], 'd'))  # (v0,v1)
        for zvals in _itertools.product(*([(0, 1)] * nqubits)):
            testvec = _functools.reduce(_np.kron, [v[i] for i in zvals])
            if _np.allclose(testvec, purevec.flat):
                return cls(zvals)
        raise ValueError(("Given `purevec` must be a z-basis product state - "
                          "cannot construct StabilizerEffectVec"))

    def __init__(self, outcomes):
        """
        Initialize a StabilizerEffectVec object.

        Parameters
        ----------
        outcomes : iterable
            A list or other iterable of integer 0 or 1 outcomes specifying
            which POVM effect vector this object represents within the
            full `stabilizerPOVM`
        """
        self._outcomes = _np.ascontiguousarray(_np.array(outcomes, int), _np.int64)
        #Note: dtype='i' => int in Cython, whereas dtype=int/np.int64 => long in Cython
        rep = replib.SBEffectRep(self._outcomes)  # dim == 2**nqubits == 2**len(outcomes)
        SPAMVec.__init__(self, rep, "stabilizer", "effect")

    #def torep(self, typ, outvec=None):
    #    # changes to_statevec/to_dmvec -> todense, and have
    #    # torep create an effect rep object...
    #    return replib.SBEffectRep(_np.ascontiguousarray(self._outcomes, _np.int64))

    def todense(self):
        """
        Return this SPAM vector as a dense state vector of shape
        (2^(nqubits), 1)

        Returns
        -------
        numpy array
        """
        v = (_np.array([1, 0], 'd'), _np.array([0, 1], 'd'))  # (v0,v1) - eigenstates of sigma_z
        statevec = _functools.reduce(_np.kron, [v[i] for i in self.outcomes])
        statevec.shape = (statevec.size, 1)
        return statevec

    @property
    def outcomes(self):
        """ The 0/1 outcomes identifying this effect within its StabilizerZPOVM """
        return self._outcomes

    def __str__(self):
        nQubits = len(self.outcomes)
        s = "Stabilizer effect vector for %d qubits with outcome %s" % (nQubits, str(self.outcomes))
        return s


class ComputationalSPAMVec(SPAMVec):
    """
    A static SPAM vector that is tensor product of 1-qubit Z-eigenstates.

    This is called a "computational basis state" in many contexts.
    """

    @classmethod
    def from_dense_vec(cls, vec, evotype):
        """
        Create a new ComputationalSPAMVec from a dense vector.

        Parameters
        ----------
        vec : numpy.ndarray
            A state vector specifying a computational basis state in the
            standard basis.  This vector has length 2^n or 4^n for
            n qubits depending on `evotype`.

        evotype : {"densitymx", "statevec", "stabilizer", "svterm", "cterm"}
            The evolution type of the resulting SPAM vector.  This value
            must be consistent with `len(vec)`, in that `"statevec"` and
            `"stabilizer"` expect 2^n whereas the rest expect 4^n.

        Returns
        -------
        ComputationalSPAMVec
        """
        if evotype in ('stabilizer', 'statevec'):
            nqubits = int(round(_np.log2(len(vec))))
            v0 = _np.array((1, 0), complex)  # '0' qubit state as complex state vec
            v1 = _np.array((0, 1), complex)  # '1' qubit state as complex state vec
        else:
            nqubits = int(round(_np.log2(len(vec)) / 2))
            v0 = 1.0 / _np.sqrt(2) * _np.array((1, 0, 0, 1), 'd')  # '0' qubit state as Pauli dmvec
            v1 = 1.0 / _np.sqrt(2) * _np.array((1, 0, 0, -1), 'd')  # '1' qubit state as Pauli dmvec

        v = (v0, v1)
        for zvals in _itertools.product(*([(0, 1)] * nqubits)):
            testvec = _functools.reduce(_np.kron, [v[i] for i in zvals])
            if _np.allclose(testvec, vec.flat):
                return cls(zvals, evotype)
        raise ValueError(("Given `vec` is not a z-basis product state - "
                          "cannot construct ComputatinoalSPAMVec"))

    def __init__(self, zvals, evotype, typ="prep"):
        """
        Initialize a ComputationalSPAMVec object.

        Parameters
        ----------
        zvals : iterable
            A list or other iterable of integer 0 or 1 outcomes specifying
            which computational basis element this object represents.  The
            length of `zvals` gives the total number of qubits.

        evotype : {"densitymx", "statevec", "stabilizer", "svterm", "cterm"}
            The type of evolution being performed.

        typ : {"prep","effect"}
            Whether this is a state preparation or POVM effect vector.
        """
        self._zvals = _np.ascontiguousarray(_np.array(zvals, _np.int64))

        nqubits = len(self._zvals)
        if evotype in ("densitymx", "svterm", "cterm"):
            dim = 4**nqubits
        elif evotype in ("statevec", "stabilizer"):
            dim = 2**nqubits
        else: raise ValueError("Invalid `evotype`: %s" % evotype)
        self._evotype = evotype  # set this before call to SPAMVec.__init__ so self.todense() can work...

        if typ == "prep":
            if evotype == "statevec":
                rep = replib.SVStateRep(self.todense())
            elif evotype == "densitymx":
                vec = _np.require(self.todense(), requirements=['OWNDATA', 'C_CONTIGUOUS'])
                rep = replib.DMStateRep(vec)
            elif evotype == "stabilizer":
                sframe = _stabilizer.StabilizerFrame.from_zvals(len(self._zvals), self._zvals)
                rep = sframe.torep()
            else:
                rep = dim  # no representations for term-based evotypes

        elif typ == "effect":
            if evotype == "statevec":
                rep = replib.SVEffectRep_Computational(self._zvals, dim)
            elif evotype == "densitymx":
                rep = replib.DMEffectRep_Computational(self._zvals, dim)
            elif evotype == "stabilizer":
                rep = replib.SBEffectRep(self._zvals)
            else:
                rep = dim  # no representations for term-based evotypes
        else:
            raise ValueError("Invalid `typ` argument for torep(): %s" % typ)

        SPAMVec.__init__(self, rep, evotype, typ)

    def todense(self, scratch=None):
        """
        Return this SPAM vector as a (dense) numpy array.  The memory
        in `scratch` maybe used when it is not-None.
        """
        if self._evotype == "densitymx":
            v0 = 1.0 / _np.sqrt(2) * _np.array((1, 0, 0, 1), 'd')  # '0' qubit state as Pauli dmvec
            v1 = 1.0 / _np.sqrt(2) * _np.array((1, 0, 0, -1), 'd')  # '1' qubit state as Pauli dmvec
        elif self._evotype in ("statevec", "stabilizer"):
            v0 = _np.array((1, 0), complex)  # '0' qubit state as complex state vec
            v1 = _np.array((0, 1), complex)  # '1' qubit state as complex state vec
        elif self._evotype in ("svterm", "cterm"):
            raise NotImplementedError("todense() is not implemented for evotype %s!" %
                                      self._evotype)
        else: raise ValueError("Invalid `evotype`: %s" % self._evotype)

        v = (v0, v1)
        return _functools.reduce(_np.kron, [v[i] for i in self._zvals])

    #def torep(self, typ, outvec=None):
    #    if typ == "prep":
    #        if self._evotype == "statevec":
    #            return replib.SVStateRep(self.todense())
    #        elif self._evotype == "densitymx":
    #            return replib.DMStateRep(self.todense())
    #        elif self._evotype == "stabilizer":
    #            sframe = _stabilizer.StabilizerFrame.from_zvals(len(self._zvals), self._zvals)
    #            return sframe.torep()
    #        raise NotImplementedError("torep(%s) not implemented for %s objects!" %
    #                                  (self._evotype, self.__class__.__name__))
    #    elif typ == "effect":
    #        if self._evotype == "statevec":
    #            return replib.SVEffectRep_Computational(self._zvals, self.dim)
    #        elif self._evotype == "densitymx":
    #            return replib.DMEffectRep_Computational(self._zvals, self.dim)
    #        elif self._evotype == "stabilizer":
    #            return replib.SBEffectRep(_np.ascontiguousarray(self._zvals, _np.int64))
    #        raise NotImplementedError("torep(%s) not implemented for %s objects!" %
    #                                  (self._evotype, self.__class__.__name__))
    #    else:
    #        raise ValueError("Invalid `typ` argument for torep(): %s" % typ)

    def get_taylor_order_terms(self, order, max_poly_vars=100, return_coeff_polys=False):
        """
        Get the `order`-th order Taylor-expansion terms of this SPAM vector.

        This function either constructs or returns a cached list of the terms at
        the given order.  Each term is "rank-1", meaning that it is a state
        preparation followed by or POVM effect preceded by actions on a
        density matrix `rho` of the form:

        `rho -> A rho B`

        The coefficients of these terms are typically polynomials of the
        SPAMVec's parameters, where the polynomial's variable indices index the
        *global* parameters of the SPAMVec's parent (usually a :class:`Model`)
        , not the SPAMVec's local parameter array (i.e. that returned from
        `to_vector`).


        Parameters
        ----------
        order : int
            The order of terms to get.

        return_coeff_polys : bool
            Whether a parallel list of locally-indexed (using variable indices
            corresponding to *this* object's parameters rather than its parent's)
            polynomial coefficients should be returned as well.

        Returns
        -------
        terms : list
            A list of :class:`RankOneTerm` objects.

        coefficients : list
            Only present when `return_coeff_polys == True`.
            A list of *compact* polynomial objects, meaning that each element
            is a `(vtape,ctape)` 2-tuple formed by concatenating together the
            output of :method:`Polynomial.compact`.
        """
        if order == 0:  # only 0-th order term exists
            if self._evotype == "svterm":
                purevec = ComputationalSPAMVec(self._zvals, "statevec", self._prep_or_effect)
            elif self._evotype == "cterm":
                purevec = ComputationalSPAMVec(self._zvals, "stabilizer", self._prep_or_effect)
            else: raise ValueError("Invalid evolution type %s for calling `get_taylor_order_terms`" % self._evotype)

            coeff = _Polynomial({(): 1.0}, max_poly_vars)
            if self._prep_or_effect == "prep":
                terms = [_term.RankOnePolyPrepTerm.simple_init(coeff, purevec, purevec, self._evotype)]
            else:
                terms = [_term.RankOnePolyEffectTerm.simple_init(coeff, purevec, purevec, self._evotype)]

            if return_coeff_polys:
                coeffs_as_compact_polys = coeff.compact(complex_coeff_tape=True)
                return terms, coeffs_as_compact_polys
            else:
                return terms  # Cache terms in FUTURE?
        else:
            if return_coeff_polys:
                vtape = _np.empty(0, _np.int64)
                ctape = _np.empty(0, complex)
                return [], (vtape, ctape)
            else:
                return []

    def num_params(self):
        """
        Get the number of independent parameters which specify this SPAM vector.

        Returns
        -------
        int
           the number of independent parameters.
        """
        return 0  # no parameters

    def to_vector(self):
        """
        Get the SPAM vector parameters as an array of values.

        Returns
        -------
        numpy array
            The parameters as a 1D array with length num_params().
        """
        return _np.array([], 'd')  # no parameters

    def from_vector(self, v, close=False, nodirty=False):
        """
        Initialize the SPAM vector using a 1D array of parameters.

        Parameters
        ----------
        v : numpy array
            The 1D vector of gate parameters.  Length
            must == num_params()

        Returns
        -------
        None
        """
        assert(len(v) == 0)  # should be no parameters, and nothing to do

    def __str__(self):
        nQubits = len(self._zvals)
        s = "Computational Z-basis SPAM vec for %d qubits w/z-values: %s" % (nQubits, str(self._zvals))
        return s