""" Utility functions operating on operation matrices """
#***************************************************************************************************
# 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.linalg as _spl
import scipy.sparse as _sps
import scipy.sparse.linalg as _spsl
import warnings as _warnings
import collections as _collections

from . import jamiolkowski as _jam
from . import matrixtools as _mt
from . import lindbladtools as _lt
from . import basistools as _bt
from ..objects.basis import Basis as _Basis, ExplicitBasis as _ExplicitBasis, DirectSumBasis as _DirectSumBasis
from ..objects.label import Label as _Label


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


def _flat_mut_blks(i, j, blockDims):
    # like _mut(i,j,dim).flatten() but works with basis *blocks*
    N = sum(blockDims)
    mx = _np.zeros((N, N), 'd'); mx[i, j] = 1.0
    ret = _np.zeros(sum([d**2 for d in blockDims]), 'd')
    i = 0; off = 0
    for d in blockDims:
        ret[i:i + d**2] = mx[off:off + d, off:off + d].flatten()
        i += d**2; off += d
    return ret


def _hack_sqrtm(A):
    sqrt, _ = _spl.sqrtm(A, disp=False)  # Travis found this scipy function
    # to be incorrect in certain cases (we need a workaround)
    if _np.any(_np.isnan(sqrt)):  # this is sometimes a good fallback when sqrtm doesn't work.
        ev, U = _np.linalg.eig(A)
        sqrt = _np.dot(U, _np.dot(_np.diag(_np.sqrt(ev)), _np.linalg.inv(U)))

    return sqrt


def fidelity(A, B):
    """
    Returns the quantum state fidelity between density
      matrices A and B given by :

      F = Tr( sqrt{ sqrt(A) * B * sqrt(A) } )^2

    To compute process fidelity, pass this function the
    Choi matrices of the two processes, or just call
    :function:`entanglement_fidelity` with the operation matrices.

    Parameters
    ----------
    A : numpy array
        First density matrix.

    B : numpy array
        Second density matrix.

    Returns
    -------
    float
        The resulting fidelity.
    """
    evals, U = _np.linalg.eig(A)
    if len([ev for ev in evals if abs(ev) > 1e-8]) == 1:
        # special case when A is rank 1, A = vec * vec^T and sqrt(A) = A
        ivec = _np.argmax(evals)
        vec = U[:, ivec:(ivec + 1)]
        F = evals[ivec].real * _np.dot(_np.conjugate(_np.transpose(vec)), _np.dot(B, vec)).real  # vec^T * B * vec
        return float(F)

    evals, U = _np.linalg.eig(B)
    if len([ev for ev in evals if abs(ev) > 1e-8]) == 1:
        # special case when B is rank 1 (recally fidelity is sym in args)
        ivec = _np.argmax(evals)
        vec = U[:, ivec:(ivec + 1)]
        F = evals[ivec].real * _np.dot(_np.conjugate(_np.transpose(vec)), _np.dot(A, vec)).real  # vec^T * A * vec
        return float(F)

    #if _np.array_equal(A, B): return 1.0  # HACK - some cases when A and B are perfecty equal sqrtm(A) fails...
    sqrtA = _hack_sqrtm(A)  # _spl.sqrtm(A)
    # test the scipy sqrtm function - sometimes fails when rank defficient
    #assert(_np.linalg.norm(_np.dot(sqrtA, sqrtA) - A) < 1e-8)
    if _np.linalg.norm(_np.dot(sqrtA, sqrtA) - A) > 1e-8:
        evals = _np.linalg.eigvals(A)
        _warnings.warn(("sqrtm(A) failure when computing fidelity - beware result. "
                        "Maybe due to rank defficiency - eigenvalues of A are: %s") % evals)
    F = (_mt.trace(_hack_sqrtm(_np.dot(sqrtA, _np.dot(B, sqrtA)))).real)**2  # Tr( sqrt{ sqrt(A) * B * sqrt(A) } )^2
    return float(F)


def frobeniusdist(A, B):
    """
    Returns the frobenius distance between gate
      or density matrices A and B given by :

      sqrt( sum( (A_ij-B_ij)^2 ) )

    Parameters
    ----------
    A : numpy array
        First matrix.

    B : numpy array
        Second matrix.

    Returns
    -------
    float
        The resulting frobenius distance.
    """
    return _mt.frobeniusnorm(A - B)


def frobeniusdist2(A, B):
    """
    Returns the square of the frobenius distance between gate
      or density matrices A and B given by :

      sum( (A_ij-B_ij)^2 )

    Parameters
    ----------
    A : numpy array
        First matrix.

    B : numpy array
        Second matrix.

    Returns
    -------
    float
        The resulting frobenius distance.
    """
    return _mt.frobeniusnorm2(A - B)


def residuals(A, B):
    """
    Calculate residuals between the elements of two matrices

    Parameters
    ----------
    A : numpy array
        First matrix.

    B : numpy array
        Second matrix.

    Returns
    -------
    np.array
        residuals
    """
    return (A - B).flatten()


def tracenorm(A):
    """
    Compute the trace norm of matrix A given by:

      Tr( sqrt{ A^dagger * A } )

    Parameters
    ----------
    A : numpy array
        The matrix to compute the trace norm of.
    """
    if _np.linalg.norm(A - _np.conjugate(A.T)) < 1e-8:
        #Hermitian, so just sum eigenvalue magnitudes
        return _np.sum(_np.abs(_np.linalg.eigvals(A)))
    else:
        #Sum of singular values (positive by construction)
        return _np.sum(_np.linalg.svd(A, compute_uv=False))


def tracedist(A, B):
    """
    Compute the trace distance between matrices A and B,
    given by:

      D = 0.5 * Tr( sqrt{ (A-B)^dagger * (A-B) } )

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the distance between.
    """
    return 0.5 * tracenorm(A - B)


def diamonddist(A, B, mxBasis='pp', return_x=False):
    """
    Returns the approximate diamond norm describing the difference between gate
    matrices A and B given by :

      D = ||A - B ||_diamond = sup_rho || AxI(rho) - BxI(rho) ||_1

    Parameters
    ----------
    A, B : numpy array
        The *gate* matrices to use when computing the diamond norm.

    mxBasis : 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).

    return_x : bool, optional
        Whether to return a numpy array encoding the state (rho) at
        which the maximal trace distance occurs.

    Returns
    -------
    dm : float
       Diamond norm
    W : numpy array
       Only returned if `return_x = True`.  Encodes the state rho, such that
       `dm = trace( |(J(A)-J(B)).T * W| )`.
    """
    mxBasis = _bt.build_basis_for_matrix(A, mxBasis)

    #currently cvxpy is only needed for this function, so don't import until here

    import cvxpy as _cvxpy

    #Check if using version < 1.0
    old_cvxpy = bool(tuple(map(int, _cvxpy.__version__.split('.'))) < (1, 0))

    # This SDP implementation is a modified version of Kevin's code

    #Compute the diamond norm

    #Uses the primal SDP from arXiv:1207.5726v2, Sec 3.2

    #Maximize 1/2 ( < J(phi), X > + < J(phi).dag, X.dag > )
    #Subject to  [[ I otimes rho0, X],
    #            [X.dag, I otimes rho1]] >> 0
    #              rho0, rho1 are density matrices
    #              X is linear operator

    #Jamiolkowski representation of the process
    #  J(phi) = sum_ij Phi(Eij) otimes Eij

    #< A, B > = Tr(A.dag B)

    #def vec(matrix_in):
    #    # Stack the columns of a matrix to return a vector
    #    return _np.transpose(matrix_in).flatten()
    #
    #def unvec(vector_in):
    #    # Slice a vector into columns of a matrix
    #    d = int(_np.sqrt(vector_in.size))
    #    return _np.transpose(vector_in.reshape( (d,d) ))

    #Code below assumes *un-normalized* Jamiol-isomorphism, so multiply by
    # density mx dimension (`smallDim`) below
    JAstd = _jam.fast_jamiolkowski_iso_std(A, mxBasis)
    JBstd = _jam.fast_jamiolkowski_iso_std(B, mxBasis)

    #Do this *after* the fast_jamiolkowski_iso calls above because these will convert
    # A & B to a "single-block" basis representation when mxBasis has multiple blocks.
    dim = JAstd.shape[0]
    smallDim = int(_np.sqrt(dim))
    JAstd *= smallDim  # see above comment
    JBstd *= smallDim  # see above comment
    assert(dim == JAstd.shape[1] == JBstd.shape[0] == JBstd.shape[1])

    #CHECK: Kevin's jamiolowski, which implements the un-normalized isomorphism:
    #  smallDim * _jam.jamiolkowski_iso(M, "std", "std")
    #def kevins_jamiolkowski(process, representation = 'superoperator'):
    #    # Return the Choi-Jamiolkowski representation of a quantum process
    #    # Add methods as necessary to accept different representations
    #    process = _np.array(process)
    #    if representation == 'superoperator':
    #        # Superoperator is the linear operator acting on vec(rho)
    #        dimension = int(_np.sqrt(process.shape[0]))
    #        print "dim = ",dimension
    #        jamiolkowski_matrix = _np.zeros([dimension**2, dimension**2], dtype='complex')
    #        for i in range(dimension**2):
    #            Ei_vec= _np.zeros(dimension**2)
    #            Ei_vec[i] = 1
    #            output = unvec(_np.dot(process,Ei_vec))
    #            tmp = _np.kron(output, unvec(Ei_vec))
    #            print "E%d = \n" % i,unvec(Ei_vec)
    #            #print "contrib =",_np.kron(output, unvec(Ei_vec))
    #            jamiolkowski_matrix += tmp
    #        return jamiolkowski_matrix
    #JAstd_kev = jamiolkowski(A)
    #JBstd_kev = jamiolkowski(B)
    #print "diff A = ",_np.linalg.norm(JAstd_kev/2.0-JAstd)
    #print "diff B = ",_np.linalg.norm(JBstd_kev/2.0-JBstd)

    #Kevin's function: def diamondnorm( jamiolkowski_matrix ):
    jamiolkowski_matrix = JBstd - JAstd

    # Here we define a bunch of auxiliary matrices because CVXPY doesn't use complex numbers

    K = jamiolkowski_matrix.real  # J.real
    L = jamiolkowski_matrix.imag  # J.imag

    if old_cvxpy:
        Y = _cvxpy.Variable(dim, dim)  # X.real
        Z = _cvxpy.Variable(dim, dim)  # X.imag

        sig0 = _cvxpy.Variable(smallDim, smallDim)  # rho0.real
        sig1 = _cvxpy.Variable(smallDim, smallDim)  # rho1.real
        tau0 = _cvxpy.Variable(smallDim, smallDim)  # rho1.imag
        tau1 = _cvxpy.Variable(smallDim, smallDim)  # rho1.imag

    else:
        Y = _cvxpy.Variable(shape=(dim, dim))  # X.real
        Z = _cvxpy.Variable(shape=(dim, dim))  # X.imag

        sig0 = _cvxpy.Variable(shape=(smallDim, smallDim))  # rho0.real
        sig1 = _cvxpy.Variable(shape=(smallDim, smallDim))  # rho1.real
        tau0 = _cvxpy.Variable(shape=(smallDim, smallDim))  # rho1.imag
        tau1 = _cvxpy.Variable(shape=(smallDim, smallDim))  # rho1.imag

    ident = _np.identity(smallDim, 'd')

    objective = _cvxpy.Maximize(_cvxpy.trace(K.T * Y + L.T * Z))
    constraints = [_cvxpy.bmat([
        [_cvxpy.kron(ident, sig0), Y, -_cvxpy.kron(ident, tau0), -Z],
        [Y.T, _cvxpy.kron(ident, sig1), Z.T, -_cvxpy.kron(ident, tau1)],
        [_cvxpy.kron(ident, tau0), Z, _cvxpy.kron(ident, sig0), Y],
        [-Z.T, _cvxpy.kron(ident, tau1), Y.T, _cvxpy.kron(ident, sig1)]]) >> 0,
        _cvxpy.bmat([[sig0, -tau0],
                     [tau0, sig0]]) >> 0,
        _cvxpy.bmat([[sig1, -tau1],
                     [tau1, sig1]]) >> 0,
        sig0 == sig0.T,
        sig1 == sig1.T,
        tau0 == -tau0.T,
        tau1 == -tau1.T,
        _cvxpy.trace(sig0) == 1.,
        _cvxpy.trace(sig1) == 1.]

    prob = _cvxpy.Problem(objective, constraints)
    try:
        prob.solve(solver="CVXOPT")
#       prob.solve(solver="ECOS")
#       prob.solve(solver="SCS")#This always fails
    except _cvxpy.error.SolverError as e:
        _warnings.warn("CVXPY failed: %s - diamonddist returning -2!" % str(e))
        return (-2, _np.zeros((dim, dim))) if return_x else -2
    except:
        _warnings.warn("CVXOPT failed (uknown err) - diamonddist returning -2!")
        return (-2, _np.zeros((dim, dim))) if return_x else -2

    #Validate result
    #assert( abs(_np.trace(_np.dot(K.T,Y.value) + _np.dot(L.T,Z.value))-prob.value) < 1e-6 ), \
    #    "Diamondnorm mismatch"

    if return_x:
        X = Y.value + 1j * Z.value  # encodes state at which maximum trace-distance occurs
        return prob.value, X
    else:
        return prob.value


def jtracedist(A, B, mxBasis='pp'):  # Jamiolkowski trace distance:  Tr(|J(A)-J(B)|)
    """
    Compute the Jamiolkowski trace distance between operation matrices A and B,
    given by:

      D = 0.5 * Tr( sqrt{ (J(A)-J(B))^2 } )

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the distance between.

    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).
    """
    JA = _jam.fast_jamiolkowski_iso_std(A, mxBasis)
    JB = _jam.fast_jamiolkowski_iso_std(B, mxBasis)
    return tracedist(JA, JB)


def entanglement_fidelity(A, B, mxBasis='pp'):
    """
    Returns the "entanglement" process fidelity between gate
    matrices A and B given by :

      F = Tr( sqrt{ sqrt(J(A)) * J(B) * sqrt(J(A)) } )^2

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the fidelity between.

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The basis of the matrices.  Allowed values are Matrix-unit (std),
        Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt)
        (or a custom basis object).
    """
    d2 = A.shape[0]
    def isTP(x): return _np.isclose(x[0, 0], 1.0) and all(
        [_np.isclose(x[0, i], 0) for i in range(d2)])

    def isUnitary(x): return _np.allclose(_np.identity(d2, 'd'), _np.dot(x, x.conjugate().T))

    if isTP(A) and isTP(B) and isUnitary(B):  # then assume TP-like gates & use simpler formula
        TrLambda = _np.trace(_np.dot(A, B.conjugate().T))  # same as using _np.linalg.inv(B)
        d2 = A.shape[0]
        return TrLambda / d2

    JA = _jam.jamiolkowski_iso(A, mxBasis, mxBasis)
    JB = _jam.jamiolkowski_iso(B, mxBasis, mxBasis)
    return fidelity(JA, JB)


def average_gate_fidelity(A, B, mxBasis='pp'):
    """
    Computes the average gate fidelity (AGF) between two gates.
    Average gate fidelity (F_g) is related to entanglement fidelity
    (F_p), via:

      F_g = (d * F_p + 1)/(1 + d),

    where d is the Hilbert space dimension. This formula, and the
    definition of AGF, can be found in Phys. Lett. A 303 249-252 (2002).

    Parameters
    ----------
    A : array or gate
        The gate to compute the AGI to B of. E.g., an imperfect
        implementation of B.

    B : array or gate
        The gate to compute the AGI to A of. E.g., the target gate
        corresponding to A.

    mxBasis : {"std","gm","pp"} or Basis object, optional
        The basis of the matrices.

    Returns
    -------
    AGI : float
        The AGI of A to B.
    """
    d = int(round(_np.sqrt(A.shape[0])))
    PF = entanglement_fidelity(A, B, mxBasis=mxBasis)
    AGF = (d * PF + 1) / (1 + d)
    return float(AGF)


def average_gate_infidelity(A, B, mxBasis="gm"):
    """
    Computes the average gate infidelity (AGI) between two gates.
    Average gate infidelity is related to entanglement infidelity
    (EI) via:

      AGI = (d * (1-EI) + 1)/(1 + d),

    where d is the Hilbert space dimension. This formula, and the
    definition of AGI, can be found in Phys. Lett. A 303 249-252 (2002).

    Parameters
    ----------
    A : array or gate
        The gate to compute the AGI to B of. E.g., an imperfect
        implementation of B.

    B : array or gate
        The gate to compute the AGI to A of. E.g., the target gate
        corresponding to A.

    mxBasis : {"std","gm","pp"} or Basis object, optional
        The basis of the matrices.

    Returns
    ----------
    AGI : float
        The AGI of A to B.
    """
    return 1 - average_gate_fidelity(A, B, mxBasis)


def entanglement_infidelity(A, B, mxBasis='pp'):
    """
    Returns the entanglement infidelity (EI) between gate
    matrices A and B given by :

      EI = 1 - Tr( sqrt{ sqrt(J(A)) * J(B) * sqrt(J(A)) } )^2

    where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix
    to it's corresponding Choi Matrix.

    Parameters
    ----------
    A, B : numpy array
        The matrices to compute the fidelity between.

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

    Returns
    -------
    EI : float
        The EI of A to B.
    """
    return 1 - float(entanglement_fidelity(A, B, mxBasis))


def gateset_infidelity(mdl, target_model, itype='EI',
                       weights=None, mxBasis=None):
    """
    Computes the average-over-gates of the infidelity between  gates in `mdl`
    and the gates in `target_model`. If `itype` is 'EI' then the "infidelity"
    is the entanglement infidelity; if `itype` is 'AGI' then the "infidelity"
    is the average gate infidelity (AGI and EI are related by a dimension
    dependent constant).

    This is the quantity that RB error rates are sometimes claimed to be
    related to directly related.

    Parameters
    ----------
    mdl : Model
        The model to calculate the average infidelity, to `target_model`, of.

    target_model : Model
        The model to calculate the average infidelity, to `mdl`, of.

    itype : str, optional
        The infidelity type. Either 'EI', corresponding to entanglement
        infidelity, or 'AGI', corresponding to average gate infidelity.

    weights : dict, optional
        If not None, a dictionary of floats, whereby the keys are the gates
        in `mdl` and the values are, possibly unnormalized, probabilities.
        These probabilities corresponding to the weighting in the average,
        so if the model contains gates A and B and weights[A] = 2 and
        weights[B] = 1 then the output is Inf(A)*2/3  + Inf(B)/3 where
        Inf(X) is the infidelity (to the corresponding element in the other
        model) of X. If None, a uniform-average is taken, equivalent to
        setting all the weights to 1.

    mxBasis : {"std","gm","pp"} or Basis object, optional
        The basis of the models. If None, the basis is obtained from
        the model.

    Returns
    -------
    float
        The weighted average-over-gates infidelity between the two models.

    """
    assert(itype == 'AGI' or itype == 'EI'), \
        "The infidelity type must be `AGI` (average gate infidelity) or `EI` (entanglement infidelity)"

    if mxBasis is None: mxBasis = mdl.basis

    sum_of_weights = 0
    I_list = []
    for gate in list(target_model.operations.keys()):
        if itype == 'AGI':
            I = average_gate_infidelity(mdl.operations[gate], target_model.operations[gate], mxBasis=mxBasis)
        if itype == 'EI':
            I = entanglement_infidelity(mdl.operations[gate], target_model.operations[gate], mxBasis=mxBasis)
        if weights is None:
            w = 1
        else:
            w = weights[gate]

        I_list.append(w * I)
        sum_of_weights += w

    assert(sum_of_weights > 0), "The sum of the weights should be positive!"
    AI = _np.sum(I_list) / sum_of_weights

    return AI


def unitarity(A, mxBasis="gm"):
    """
    Returns the "unitarity" of a channel, as defined in Wallman et al,
    ``Estimating the Coherence of noise'' NJP 17 113020 (2015). The
    unitarity is given by (Prop 1 in Wallman et al):

    u(A) = Tr( A_u^{\dagger} A_u ) / (d^2  - 1),

    where A_u is the unital submatrix of A, and d is the dimension of
    the Hilbert space. When A is written in any basis for which the
    first element is the  normalized identity (e.g., the pp or gm
    bases), The unital submatrix of A is the matrix obtained when the
    top row and left hand column is removed from A.

    Parameters
    ----------
    A : array or gate
        The gate for which the unitarity is to be computed.

    mxBasis : {"std","gm","pp"} or a Basis object, optional
        The basis of the matrix.

    d : int, optional
        The dimension of the Hilbert space.

    Returns
    ----------
    u : float
        The unitarity of the gate A.

    """
    d = int(round(_np.sqrt(A.shape[0])))
    basisMxs = _bt.basis_matrices(mxBasis, A.shape[0])

    if _np.allclose(basisMxs[0], _np.identity(d, 'd')):
        B = A
    else:
        B = _bt.change_basis(A, mxBasis, "gm")  # everything should be able to be put in the "gm" basis

    unital = B[1:d**2, 1:d**2]
    u = _np.trace(_np.dot(_np.conj(_np.transpose(unital)), unital)) / (d**2 - 1)
    return u


def fidelity_upper_bound(operationMx):
    """
    Get an upper bound on the fidelity of the given
      operation matrix with any unitary operation matrix.

    The closeness of the result to one tells
     how "unitary" the action of operationMx is.

    Parameters
    ----------
    operationMx : numpy array
        The operation matrix to act on.

    Returns
    -------
    float
        The resulting upper bound on fidelity(operationMx, anyUnitaryGateMx)
    """
    choi = _jam.jamiolkowski_iso(operationMx, choiMxBasis="std")
    choi_evals, choi_evecs = _np.linalg.eig(choi)
    maxF_direct = max([_np.sqrt(max(ev.real, 0.0)) for ev in choi_evals]) ** 2

    iMax = _np.argmax([ev.real for ev in choi_evals])  # index of maximum eigenval
    closestVec = choi_evecs[:, iMax:(iMax + 1)]

    # #print "DEBUG: closest evec = ", closestUnitaryVec
    # new_evals = _np.zeros( len(closestUnitaryVec) ); new_evals[iClosestU] = 1.0
    # # gives same result:
    # closestUnitaryJmx = _np.dot(choi_evecs, _np.dot( _np.diag(new_evals), _np.linalg.inv(choi_evecs) ) )
    closestJmx = _np.kron(closestVec, _np.transpose(_np.conjugate(closestVec)))  # closest rank-1 Jmx
    closestJmx /= _mt.trace(closestJmx)  # normalize so trace of Jmx == 1.0

    maxF = fidelity(choi, closestJmx)

    if not _np.isnan(maxF):

        #Uncomment for debugging
        #if abs(maxF - maxF_direct) >= 1e-6:
        #    print "DEBUG: operationMx:\n",operationMx
        #    print "DEBUG: choiMx:\n",choi
        #    print "DEBUG choi_evals = ",choi_evals, " iMax = ",iMax
        #    #print "DEBUG: J = \n", closestUnitaryJmx
        #    print "DEBUG: eigvals(J) = ", _np.linalg.eigvals(closestJmx)
        #    print "DEBUG: trace(J) = ", _mt.trace(closestJmx)
        #    print "DEBUG: maxF = %f,  maxF_direct = %f" % (maxF, maxF_direct)
        #    raise ValueError("ERROR: maxF - maxF_direct = %f" % (maxF -maxF_direct))
        assert(abs(maxF - maxF_direct) < 1e-6)
    else:
        maxF = maxF_direct  # case when maxF is nan, due to scipy sqrtm function being buggy - just use direct F

    closestOpMx = _jam.jamiolkowski_iso_inv(closestJmx, choiMxBasis="std")
    return maxF, closestOpMx

    #closestU_evals, closestU_evecs = _np.linalg.eig(closestUnitaryGateMx)
    #print "DEBUG: U = \n", closestUnitaryGateMx
    #print "DEBUG: closest U evals = ",closestU_evals
    #print "DEBUG:  evecs = \n",closestU_evecs


def get_povm_map(model, povmlbl):
    """
    Constructs a gate-like quantity for the POVM within `model`.

    This is done by embedding the `k`-outcome classical output space of the POVM
    in the Hilbert-Schmidt space of `k` by `k` density matrices by placing the
    classical probability distribution along the diagonal of the density matrix.
    Currently, this is only implemented for the case when `k` equals `d`, the
    dimension of the POVM's Hilbert space.

    Parameters
    ----------
    model : Model
        The model supplying the POVM effect vectors and the basis those
        vectors are in.

    povmlbl : str
        The POVM label

    Returns
    -------
    numpy.ndarray
        The matrix of the "POVM map" in the `model.basis` basis.
    """
    povmVectors = [v.todense()[:, None] for v in model.povms[povmlbl].values()]
    if isinstance(model.basis, _DirectSumBasis):  # HACK - need to get this to work with general bases
        blkDims = [int(_np.sqrt(comp.dim)) for comp in model.basis.component_bases]
    else:
        blkDims = [int(round(_np.sqrt(model.dim)))]  # [d] where density matrix is dxd

    nV = len(povmVectors)
    #assert(d**2 == model.dim), "Model dimension (%d) is not a perfect square!" % model.dim
    #assert( nV**2 == d ), "Can only compute POVM metrics when num of effects == H space dimension"
    #   I don't think above assert is needed - should work in general (Robin?)
    povm_mx = _np.concatenate(povmVectors, axis=1).T  # "povm map" ( B(H) -> S_k ) (shape= nV,model.dim)

    Sk_embedding_in_std = _np.zeros((model.dim, nV))
    for i in range(nV):
        Sk_embedding_in_std[:, i] = _flat_mut_blks(i, i, blkDims)

    std_to_basis = model.basis.reverse_transform_matrix("std")  # _bt.transform_matrix("std", model.basis, blkDims)
    assert(std_to_basis.shape == (model.dim, model.dim))

    return _np.dot(std_to_basis, _np.dot(Sk_embedding_in_std, povm_mx))


def povm_fidelity(model, targetModel, povmlbl):
    """
    Computes the process (entanglement) fidelity between POVM maps.

    Parameters
    ----------
    model, targetModel : Model
        Models containing the two POVMs to compare.

    povmlbl : str
        The POVM label

    Returns
    -------
    float
    """
    povm_mx = get_povm_map(model, povmlbl)
    target_povm_mx = get_povm_map(targetModel, povmlbl)
    return entanglement_fidelity(povm_mx, target_povm_mx, targetModel.basis)


def povm_jtracedist(model, targetModel, povmlbl):
    """
    Computes the Jamiolkowski trace distance between POVM maps using :func:`jtracedist`.

    Parameters
    ----------
    model, targetModel : Model
        Models containing the two POVMs to compare.

    povmlbl : str
        The POVM label

    Returns
    -------
    float
    """
    povm_mx = get_povm_map(model, povmlbl)
    target_povm_mx = get_povm_map(targetModel, povmlbl)
    return jtracedist(povm_mx, target_povm_mx, targetModel.basis)


def povm_diamonddist(model, targetModel, povmlbl):
    """
    Computes the diamond distance between POVM maps using :func:`diamonddist`.

    Parameters
    ----------
    model, targetModel : Model
        Models containing the two POVMs to compare.

    povmlbl : str
        The POVM label

    Returns
    -------
    float
    """
    povm_mx = get_povm_map(model, povmlbl)
    target_povm_mx = get_povm_map(targetModel, povmlbl)
    return diamonddist(povm_mx, target_povm_mx, targetModel.basis)


#decompose operation matrix into axis of rotation, etc
def decompose_gate_matrix(operationMx):
    """
    Compute how the action of a operation matrix can be
    is decomposed into fixed points, axes of rotation,
    angles of rotation, and decays.  Also determines
    whether a gate appears to be valid and/or unitary.

    Parameters
    ----------
    operationMx : numpy array
        The operation matrix to act on.

    Returns
    -------
    dict
       A dictionary describing the decomposed action. Keys are:

         'isValid' : bool
             whether decomposition succeeded
         'isUnitary' : bool
             whether operationMx describes unitary action
         'fixed point' : numpy array
             the fixed point of the action
         'axis of rotation' : numpy array or nan
             the axis of rotation
         'decay of diagonal rotation terms' : float
             decay of diagonal terms
         'rotating axis 1' : numpy array or nan
             1st axis orthogonal to axis of rotation
         'rotating axis 2' : numpy array or nan
             2nd axis orthogonal to axis of rotation
         'decay of off diagonal rotation terms' : float
             decay of off-diagonal terms
         'pi rotations' : float
             angle of rotation in units of pi radians
    """

    op_evals, op_evecs = _np.linalg.eig(_np.asarray(operationMx))
    # fp_eigenvec = None
    # aor_eval = None; aor_eigenvec = None
    # ra_eval  = None; ra1_eigenvec = None; ra2_eigenvec = None

    TOL = 1e-4  # 1e-7

    unit_eval_indices = [i for (i, ev) in enumerate(op_evals) if abs(ev - 1.0) < TOL]
    #unit_eval_indices = [ i for (i,ev) in enumerate(op_evals) if ev > (1.0-TOL) ]

    conjpair_eval_indices = []
    for (i, ev) in enumerate(op_evals):
        if i in unit_eval_indices: continue  # don't include the unit eigenvalues in the conjugate pair count
        # don't include existing conjugate pairs
        if any([(i in conjpair) for conjpair in conjpair_eval_indices]): continue
        for (j, ev2) in enumerate(op_evals[i + 1:]):
            if abs(ev - _np.conjugate(ev2)) < TOL:
                conjpair_eval_indices.append((i, j + (i + 1)))
                break  # don't pair i-th eigenvalue with any other (pairs should be disjoint)

    real_eval_indices = []  # indices of real eigenvalues that are not units or a part of any conjugate pair
    complex_eval_indices = []  # indices of complex eigenvalues that are not units or a part of any conjugate pair
    for (i, ev) in enumerate(op_evals):
        if i in unit_eval_indices: continue  # don't include the unit eigenvalues
        if any([(i in conjpair) for conjpair in conjpair_eval_indices]): continue  # don't include the conjugate pairs
        if abs(ev.imag) < TOL: real_eval_indices.append(i)
        else: complex_eval_indices.append(i)

    #if len(real_eval_indices + unit_eval_indices) > 0:
    #    max_real_eval = max([ op_evals[i] for i in real_eval_indices + unit_eval_indices])
    #    min_real_eval = min([ op_evals[i] for i in real_eval_indices + unit_eval_indices])
    #else:
    #    max_real_eval = _np.nan
    #    min_real_eval = _np.nan
    #
    #fixed_points = [ op_evecs[:,i] for i in unit_eval_indices ]
    #real_eval_axes = [ op_evecs[:,i] for i in real_eval_indices ]
    #conjpair_eval_axes = [ (op_evecs[:,i],op_evecs[:,j]) for (i,j) in conjpair_eval_indices ]
    #
    #ret = { }

    nQubits = _np.log2(operationMx.shape[0]) / 2
    if nQubits == 1:
        #print "DEBUG: 1 qubit decomp --------------------------"
        #print "   --> evals = ", op_evals
        #print "   --> unit eval indices = ", unit_eval_indices
        #print "   --> conj eval indices = ", conjpair_eval_indices
        #print "   --> unpaired real eval indices = ", real_eval_indices

        #Special case: if have two conjugate pairs, check if one (or both) are real
        #  and break the one with the largest (real) value into two unpaired real evals.
        if len(conjpair_eval_indices) == 2:
            iToBreak = None
            if abs(_np.imag(op_evals[conjpair_eval_indices[0][0]])) < TOL and \
               abs(_np.imag(op_evals[conjpair_eval_indices[1][0]])) < TOL:
                iToBreak = _np.argmax([_np.real(conjpair_eval_indices[0][0]), _np.real(conjpair_eval_indices[1][0])])
            elif abs(_np.imag(op_evals[conjpair_eval_indices[0][0]])) < TOL: iToBreak = 0
            elif abs(_np.imag(op_evals[conjpair_eval_indices[1][0]])) < TOL: iToBreak = 1

            if iToBreak is not None:
                real_eval_indices.append(conjpair_eval_indices[iToBreak][0])
                real_eval_indices.append(conjpair_eval_indices[iToBreak][1])
                del conjpair_eval_indices[iToBreak]

        #Find eigenvector corresponding to fixed point (or closest we can get).   This
        # should be a unit eigenvalue with identity eigenvector.
        if len(unit_eval_indices) > 0:
            #Find linear least squares solution within possibly degenerate unit-eigenvalue eigenspace
            # of eigenvector closest to identity density mx (the desired fixed point), then orthogonalize
            # the remaining eigenvectors w.r.t this one.
            A = _np.take(op_evecs, unit_eval_indices, axis=1)
            b = _np.array([[1], [0], [0], [0]], 'd')  # identity density mx
            x = _np.dot(_np.linalg.pinv(_np.dot(A.T, A)), _np.dot(A.T, b))
            fixedPtVec = _np.dot(A, x)  # fixedPtVec / _np.linalg.norm(fixedPtVec)
            fixedPtVec = fixedPtVec[:, 0]

            iLargestContrib = _np.argmax(_np.abs(x))  # index of gate eigenvector which contributed the most
            for ii, i in enumerate(unit_eval_indices):
                if ii == iLargestContrib:
                    op_evecs[:, i] = fixedPtVec
                    iFixedPt = i
                else:
                    op_evecs[:, i] = op_evecs[:, i] - _np.vdot(fixedPtVec, op_evecs[:, i]) * fixedPtVec
                    for jj, j in enumerate(unit_eval_indices[:ii]):
                        if jj == iLargestContrib: continue
                        op_evecs[:, i] = op_evecs[:, i] - _np.vdot(op_evecs[:, j], op_evecs[:, i]) * op_evecs[:, j]
                    op_evecs[:, i] /= _np.linalg.norm(op_evecs[:, i])

        elif len(real_eval_indices) > 0:
            # just take eigenvector corresponding to the largest real eigenvalue?
            #iFixedPt = real_eval_indices[ _np.argmax( [ op_evals[i] for i in real_eval_indices ] ) ]

            # ...OR take eigenvector corresponding to a real unpaired eigenvalue closest to identity:
            idmx = _np.array([[1], [0], [0], [0]], 'd')  # identity density mx
            iFixedPt = real_eval_indices[_np.argmin([_np.linalg.norm(op_evecs[i] - idmx) for i in real_eval_indices])]

        else:
            #No unit or real eigenvalues => two complex conjugate pairs or unpaired complex evals --> bail out
            return {'isValid': False, 'isUnitary': False, 'msg': "All evals are complex."}

        #Find eigenvector corresponding to axis of rotation: find the *largest* unpaired real/unit eval
        indsToConsider = (unit_eval_indices + real_eval_indices)[:]
        del indsToConsider[indsToConsider.index(iFixedPt)]  # don't consider fixed pt evec

        if len(indsToConsider) > 0:
            iRotAxis = indsToConsider[_np.argmax([op_evals[i] for i in indsToConsider])]
        else:
            #No unit or real eigenvalues => an unpaired complex eval --> bail out
            return {'isValid': False, 'isUnitary': False, 'msg': "Unpaired complex eval."}

        #There are only 2 eigenvalues left -- hopefully a conjugate pair giving rotation
        inds = list(range(4))
        del inds[inds.index(iFixedPt)]
        del inds[inds.index(iRotAxis)]
        if abs(op_evals[inds[0]] - _np.conjugate(op_evals[inds[1]])) < TOL:
            iConjPair1, iConjPair2 = inds
        else:
            return {'isValid': False, 'isUnitary': False, 'msg': "No conjugate pair for rotn."}

        return {'isValid': True,
                'isUnitary': bool(len(unit_eval_indices) >= 2),
                'fixed point': op_evecs[:, iFixedPt],
                'axis of rotation': op_evecs[:, iRotAxis],
                'rotating axis 1': op_evecs[:, iConjPair1],
                'rotating axis 2': op_evecs[:, iConjPair2],
                'decay of diagonal rotation terms': 1.0 - abs(op_evals[iRotAxis]),
                'decay of off diagonal rotation terms': 1.0 - abs(op_evals[iConjPair1]),
                'pi rotations': _np.angle(op_evals[iConjPair1]) / _np.pi,
                'msg': "Success"}

    else:
        return {'isValid': False,
                'isUnitary': False,
                'msg': "Unsupported number of qubits: %d" % nQubits}


def state_to_dmvec(psi):
    """
    Compute the vectorized density matrix which acts as the state `psi`.

    This is just the outer product map |psi> => |psi><psi| with the
    output flattened, i.e. `dot(psi, conjugate(psi).T)`.

    Parameters
    ----------
    psi : numpy array
        The state vector.

    Returns
    -------
    numpy array
       The vectorized density matrix.
    """
    psi = psi.reshape((psi.size, 1))  # convert to (N,1) shape if necessary
    dm = _np.dot(psi, _np.conjugate(psi.T))
    return dm.flatten()


def dmvec_to_state(dmvec, tol=1e-6):
    """
    Compute the pure state describing the action of density matrix vector `dmvec`.

    If `dmvec` represents a mixed state, ValueError is raised.

    Parameters
    ----------
    dmvec : numpy array
        The vectorized density matrix, assumed to be in the standard (matrix
        unit) basis.

    tol : float, optional
        tolerance for determining whether an eigenvalue is zero.

    Returns
    -------
    numpy array
       The pure state, as a column vector of shape = (N,1)
    """
    d2 = dmvec.size; d = int(round(_np.sqrt(d2)))
    dm = dmvec.reshape((d, d))
    evals, evecs = _np.linalg.eig(dm)

    k = None
    for i, ev in enumerate(evals):
        if abs(ev) > tol:
            if k is None: k = i
            else: raise ValueError("Cannot convert mixed dmvec to pure state!")
    if k is None: raise ValueError("Cannot convert zero dmvec to puse state!")
    psi = evecs[:, k] * _np.sqrt(evals[k])
    psi.shape = (d, 1)
    return psi


def unitary_to_process_mx(U):
    """
    Compute the super-operator which acts on (row)-vectorized
    density matrices from a unitary operator (matrix) U which
    acts on state vectors.  This super-operator is given by
    the tensor product of U and conjugate(U), i.e. kron(U,U.conj).

    Parameters
    ----------
    U : numpy array
        The unitary matrix which acts on state vectors.

    Returns
    -------
    numpy array
       The super-operator process matrix.
    """
    # U -> kron(U,Uc) since U rho U_dag -> kron(U,Uc)
    #  since AXB --row-vectorize--> kron(A,B.T)*vec(X)
    return _np.kron(U, _np.conjugate(U))


def process_mx_to_unitary(superop):
    """
    Compute the unitary corresponding to the (unitary-action!)
    super-operator `superop` which acts on (row)-vectorized
    density matrices.  The super-operator must be of the form
    `kron(U,U.conj)` or an error will be thrown.

    Parameters
    ----------
    superop : numpy array
        The superoperator matrix which acts on vectorized
        density matrices (in the 'std' matrix-unit basis).

    Returns
    -------
    numpy array
       The unitary matrix which acts on state vectors.
    """
    d2 = superop.shape[0]; d = int(round(_np.sqrt(d2)))
    U = _np.empty((d, d), 'complex')

    for i in range(d):
        densitymx_i = _np.zeros((d, d), 'd'); densitymx_i[i, i] = 1.0  # |i><i|
        UiiU = _np.dot(superop, densitymx_i.flat).reshape((d, d))  # U|i><i|U^dag

        if i > 0:
            j = 0
            densitymx_ij = _np.zeros((d, d), 'd'); densitymx_ij[i, j] = 1.0  # |i><i|
            UijU = _np.dot(superop, densitymx_ij.flat).reshape((d, d))  # U|i><j|U^dag
            Uj = U[:, j]
            Ui = _np.dot(UijU, Uj)
        else:
            ##method1: use random state projection
            #rand_state = _np.random.rand(d)
            #projected_rand_state = _np.dot(UiiU, rand_state)
            #assert(_np.linalg.norm(projected_rand_state) > 1e-8)
            #projected_rand_state /= _np.linalg.norm(projected_rand_state)
            #Ui = projected_rand_state

            #method2: get eigenvector corresponding to largest eigenvalue (more robust)
            evals, evecs = _np.linalg.eig(UiiU)
            imaxeval = _np.argmax(_np.abs(evals))
            #TODO: assert other eigenvalues are much smaller?
            Ui = evecs[:, imaxeval]
            Ui /= _np.linalg.norm(Ui)
        U[:, i] = Ui

    return U


def spam_error_generator(spamvec, target_spamvec, mxBasis, typ="logGTi"):
    """
    Construct an error generator from a SPAM vector and it's target.

    Computes the value of the error generator given by
    `errgen = log( diag(spamvec / target_spamvec) )`, where division is
    element-wise.  This results in a (non-unique) error generator matrix
    `E` such that `spamvec = exp(E) * target_spamvec`.

    Note: This is currently of very limited use, as the above algorithm fails
    whenever `target_spamvec` has zero elements where `spamvec` doesn't.

    Parameters
    ----------
    spamvec : ndarray
      The SPAM vector.

    target_spamvec : ndarray
      The target SPAM vector.

    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).

    typ : {"logGTi"}
      The type of error generator to compute.  Allowed values are:

      - "logGTi" : errgen = log( diag(spamvec / target_spamvec) )

    Returns
    -------
    errgen : ndarray
      The error generator.
    """
    # Compute error generator for rho:   rho = exp(E)rho0 => rho = A*rho0 => A = diag(rho/rho0)
    assert(typ == "logGTi"), "Only logGTi type is supported so far"

    d2 = len(spamvec)
    errgen = _np.zeros((d2, d2), 'd')  # type assumes this is density-mx evolution
    diags = []
    for a, b in zip(spamvec, target_spamvec):
        if _np.isclose(b, 0.0):
            if _np.isclose(a, b): d = 1
            else: raise ValueError("Cannot take spam_error_generator")
        else:
            d = a / b
        diags.append(d)
    errgen[_np.diag_indices(d2)] = diags
    return _spl.logm(errgen)


def error_generator(gate, target_op, mxBasis, typ="logG-logT"):
    """
    Construct the error generator from a gate and its target.

    Computes the value of the error generator given by
    errgen = log( inv(target_op) * gate ), so that
    gate = target_op * exp(errgen).

    Parameters
    ----------
    gate : ndarray
      The operation matrix

    target_op : ndarray
      The target operation matrix

    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).

    typ : {"logG-logT", "logTiG", "logGTi"}
      The type of error generator to compute.  Allowed values are:

      - "logG-logT" : errgen = log(gate) - log(target_op)
      - "logTiG" : errgen = log( dot(inv(target_op), gate) )
      - "logGTi" : errgen = log( dot(gate,inv(target_op)) )

    Returns
    -------
    errgen : ndarray
      The error generator.
    """
    TOL = 1e-8

    if typ == "logG-logT":
        try:
            logT = _mt.unitary_superoperator_matrix_log(target_op, mxBasis)
        except AssertionError:  # if not unitary, fall back to just taking the real log
            logT = _mt.real_matrix_log(target_op, "raise", TOL)  # make a fuss if this can't be done
        logG = _mt.approximate_matrix_log(gate, logT)

        # Both logG and logT *should* be real, so we just take the difference.
        if _np.linalg.norm(_np.imag(logG)) < TOL and \
           _np.linalg.norm(_np.imag(logT)) < TOL:
            return _np.real(logG - logT)

        #Otherwise, there could be branch cut issues or worse, so just
        # raise an error for now (maybe return a dummy if needed elsewhere?)
        raise ValueError("Could not construct a real logarithms for the "
                         "'logG-logT' generator.  Perhaps you should use "
                         "the 'logTiG' or 'logGTi' generator instead?")

    elif typ == "logTiG":
        target_op_inv = _spl.inv(target_op)
        try:
            errgen = _mt.near_identity_matrix_log(_np.dot(target_op_inv, gate), TOL)
        except AssertionError:  # not near the identity, fall back to the real log
            _warnings.warn(("Near-identity matrix log failed; falling back "
                            "to approximate log for logTiG error generator"))
            errgen = _mt.real_matrix_log(_np.dot(target_op_inv, gate), "warn", TOL)

        if _np.linalg.norm(errgen.imag) > TOL:
            _warnings.warn("Falling back to approximate log for logTiG error generator")
            errgen = _mt.approximate_matrix_log(_np.dot(target_op_inv, gate),
                                                _np.zeros(gate.shape, 'd'), TOL=TOL)

    elif typ == "logGTi":
        target_op_inv = _spl.inv(target_op)
        try:
            errgen = _mt.near_identity_matrix_log(_np.dot(gate, target_op_inv), TOL)
        except AssertionError as e:  # not near the identity, fall back to the real log
            _warnings.warn(("Near-identity matrix log failed; falling back "
                            "to approximate log for logGTi error generator:\n%s") % str(e))
            errgen = _mt.real_matrix_log(_np.dot(gate, target_op_inv), "warn", TOL)

        if _np.linalg.norm(errgen.imag) > TOL:
            _warnings.warn("Falling back to approximate log for logGTi error generator")
            errgen = _mt.approximate_matrix_log(_np.dot(gate, target_op_inv),
                                                _np.zeros(gate.shape, 'd'), TOL=TOL)

    else:
        raise ValueError("Invalid error-generator type: %s" % typ)

    if _np.linalg.norm(_np.imag(errgen)) > TOL:
        raise ValueError("Could not construct a real generator!")
        #maybe this is actually ok, but a complex error generator will
        # need to be plotted differently, etc -- TODO
    return _np.real(errgen)


def operation_from_error_generator(error_gen, target_op, typ="logG-logT"):
    """
    Construct a gate from an error generator and a target gate.

    Inverts the computation fone in :func:`error_generator` and
    returns the value of the gate given by
    gate = target_op * exp(error_gen).

    Parameters
    ----------
    error_gen : ndarray
      The error generator matrix

    target_op : ndarray
      The target operation matrix

    typ : {"logG-logT", "logTiG"}
      The type of error generator to compute.  Allowed values are:

      - "logG-logT" : errgen = log(gate) - log(target_op)
      - "logTiG" : errgen = log( dot(inv(target_op), gate) )


    Returns
    -------
    ndarray
      The operation matrix.
    """
    if typ == "logG-logT":
        return _spl.expm(error_gen + _spl.logm(target_op))
    elif typ == "logTiG":
        return _np.dot(target_op, _spl.expm(error_gen))
    elif typ == "logGTi":
        return _np.dot(_spl.expm(error_gen), target_op)
    else:
        raise ValueError("Invalid error-generator type: %s" % typ)


def std_scale_factor(dim, projection_type):
    """
    Returns the multiplicative scaling that should be applied to the output of
    :func"`std_error_generators`, before using them as projectors, in order to
    compute the "standard" reported projection onto that type of error (i.e.
    the coefficient of the standard generator terms built un-normalized-Paulis).

    Parameters
    ----------
    dim : int
      The dimension of the error generators; also the  associated gate
      dimension.  This must be a perfect square, as `sqrt(dim)`
      is the dimension of density matrices. For a single qubit, dim == 4.

    projection_type : {"hamiltonian", "stochastic", "affine"}
      The type/class of error generators to get the scaling for.

    Returns
    -------
    float
    """
    d2 = dim
    d = int(_np.sqrt(d2))

    if projection_type == "hamiltonian":
        scaleFctr = 1.0 / (d * _np.sqrt(2))
        # so projection is coefficient of Hamiltonian term (w/un-normalized Paulis)
    elif projection_type == "stochastic":
        scaleFctr = 1.0 / d
        # so projection is coefficient of P*rho*P stochastic term in generator (w/un-normalized Paulis)
    elif projection_type == "affine":
        scaleFctr = 1.0  # so projection is coefficient of P affine term in generator (w/un-normalized Paulis)
    else:
        raise ValueError("Invalid projection_type argument: %s"
                         % projection_type)
    return scaleFctr


def std_error_generators(dim, projection_type, projection_basis):
    """
    Compute the gate error generators for a standard set of errors which
    correspond to "Hamiltonian"- or "Stochastic"-type errors in terms of the
    elements of the specified basis.

    Parameters
    ----------
    dim : int
      The dimension of the error generators to be returned.  This is also the
      associated gate dimension, and must be a perfect square, as `sqrt(dim)`
      is the dimension of density matrices. For a single qubit, dim == 4.

    projection_type : {"hamiltonian", "stochastic", "affine"}
      The type of error generators to construct.  If "hamiltonian", then the
      Hamiltonian generators which take a density matrix rho -> -i*[ H, rho ]
      for Pauli-product matrix H.  If "stochastic", then the Stochastic error
      generators which take rho -> P*rho*P for Pauli-product matrix P.  If
      "affine", then the affine generators which take rho -> P.

    projection_basis : {'std', 'gm', 'pp', 'qt'}
      Which basis is used to construct the error generators.  Allowed
      values are Matrix-unit (std), Gell-Mann (gm),
      Pauli-product (pp) and Qutrit (qt).

    Returns
    -------
    generators : numpy.ndarray
      An array of shape (#basis-elements,dim,dim).  `generators[i]` is the
      generator corresponding to the ith basis matrix in the
      *std* (matrix unit) basis.  (Note that in most cases #basis-elements
      == dim, so the size of `generators` is (dim,dim,dim) ).  Each
      generator is normalized so that as a vector it has unit Frobenius norm.
    """
    d2 = dim
    d = int(_np.sqrt(d2))

    #Get a list of the basis matrices
    mxs = _bt.basis_matrices(projection_basis, d2)

    assert(len(mxs) <= d2)  # OK if there are fewer basis matrices (e.g. for bases w/multiple blocks)
    assert(_np.isclose(d * d, d2))  # d2 must be a perfect square

    lindbladMxs = _np.empty((len(mxs), d2, d2), 'complex')
    for i, basisMx in enumerate(mxs):
        if projection_type == "hamiltonian":
            lindbladMxs[i] = _lt.hamiltonian_to_lindbladian(basisMx)  # in std basis
        elif projection_type == "stochastic":
            lindbladMxs[i] = _lt.stochastic_lindbladian(basisMx)  # in std basis
        elif projection_type == "affine":
            lindbladMxs[i] = _lt.affine_lindbladian(basisMx)  # in std basis
        else:
            raise ValueError("Invalid projection_type argument: %s"
                             % projection_type)
        norm = _np.linalg.norm(lindbladMxs[i].flat)
        if not _np.isclose(norm, 0):
            lindbladMxs[i] /= norm  # normalize projector
            assert(_np.isclose(_np.linalg.norm(lindbladMxs[i].flat), 1.0))

    return lindbladMxs


def std_errgen_projections(errgen, projection_type, projection_basis,
                           mxBasis="gm", return_generators=False,
                           return_scale_fctr=False):
    """
    Compute the projections of a gate error generator onto generators
    for a standard set of errors constructed from the elements of a
    specified basis.

    Parameters
    ----------
    errgen: : ndarray
      The error generator matrix to project.

    projection_type : {"hamiltonian", "stochastic", "affine"}
      The type of error generators to project the gate error generator onto.
      If "hamiltonian", then use the Hamiltonian generators which take a density
      matrix rho -> -i*[ H, rho ] for Pauli-product matrix H.  If "stochastic",
      then use the Stochastic error generators which take rho -> P*rho*P for
      Pauli-product matrix P (recall P is self adjoint).  If "affine", then
      use the affine error generators which take rho -> P (superop is |P>><<1|).

    projection_basis : {'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).

    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).

    return_generators : bool, optional
      If True, return the error generators projected against along with the
      projection values themseves.

    return_scale_fctr : bool, optional
      If True, also return the scaling factor that was used to multply the
      projections onto *normalized* error generators to get the returned
      values.

    Returns
    -------
    projections : numpy.ndarray
      An array of length equal to the number of elements in the
      basis used to construct the projectors.  Typically this is
      is also the dimension of the gate (e.g. 4 for a single qubit).

    generators : numpy.ndarray
      Only returned when `return_generators == True`.  An array of shape
      (#basis-els,op_dim,op_dim) such that  `generators[i]` is the
      generator corresponding to the i-th basis element.  Note
      that these matricies are in the *std* (matrix unit) basis.

    scale : float
      Only returned when `return_scale_fctr == True`.  A mulitplicative
      scaling constant that *has already been applied* to `projections`.
    """

    if isinstance(mxBasis, _Basis):
        errgen_std = _bt.change_basis(errgen, mxBasis, mxBasis.equivalent('std'))

        #expand operation matrix so it acts on entire space of dmDim x dmDim density matrices
        errgen_std = _bt.resize_std_mx(errgen_std, 'expand', mxBasis.equivalent('std'),
                                       mxBasis.simple_equivalent('std'))
    else:
        errgen_std = _bt.change_basis(errgen, mxBasis, "std")

    d2 = errgen_std.shape[0]
    d = int(_np.sqrt(d2))
    # nQubits = _np.log2(d)

    #Get a list of the d2 generators (in corresspondence with the
    #  Pauli-product matrices given by _basis.pp_matrices(d) ).
    lindbladMxs = std_error_generators(d2, projection_type, projection_basis)  # in std basis

    assert(len(lindbladMxs) <= d2)  # can be fewer projection matrices (== lenght of projection_basis)
    assert(_np.isclose(d * d, d2))  # d2 must be a perfect square

    projections = _np.empty(len(lindbladMxs), 'd')
    for i, lindbladMx in enumerate(lindbladMxs):
        proj = _np.real_if_close(_np.vdot(errgen_std.flatten(), lindbladMx.flatten()), tol=1000)

        # # DEBUG - for checking why perfect gates gave weird projections --> log ambiguity
        # print("DB: rawproj(%d) = " % i, proj)
        # errgen_pp = errgen.copy() #_bt.change_basis(errgen_std,"std","pp")
        # lindbladMx_pp = _bt.change_basis(lindbladMx,"std","pp")
        # if proj > 1.0:
        #    for k in range(errgen_std.shape[0]):
        #        for j in range(errgen_std.shape[1]):
        #            if abs(errgen_pp[k,j].conjugate() * lindbladMx_pp[k,j]) > 1e-2:
        #                print(" [%d,%d]: + " % (k,j), errgen_pp[k,j].conjugate(),
        #                      "*", lindbladMx_pp[k,j],
        #                      "=", (errgen_pp[k,j].conjugate() * lindbladMx_pp[i,j]))

        #assert(_np.isreal(proj)), "non-real projection: %s" % str(proj) #just a warning now
        if not _np.isreal(proj):
            _warnings.warn("Taking abs() of non-real projection: %s" % str(proj))
            proj = abs(proj)
        projections[i] = proj

    scaleFctr = std_scale_factor(d2, projection_type)
    projections *= scaleFctr
    lindbladMxs /= scaleFctr  # so projections * generators give original

    ret = [projections]
    if return_generators: ret.append(lindbladMxs)
    if return_scale_fctr: ret.append(scaleFctr)
    return ret[0] if len(ret) == 1 else tuple(ret)


def _assert_shape(ar, shape, sparse=False):
    """ Asserts ar.shape == shape ; works with sparse matrices too """
    if not sparse or len(shape) == 2:
        assert(ar.shape == shape), \
            "Shape mismatch: %s != %s!" % (str(ar.shape), str(shape))
    else:
        if len(shape) == 3:  # first "dim" is a list
            assert(len(ar) == shape[0]), \
                "Leading dim mismatch: %d != %d!" % (len(ar), shape[0])
            assert(shape[0] == 0 or ar[0].shape == (shape[1], shape[2])), \
                "Shape mismatch: %s != %s!" % (str(ar[0].shape), str(shape[1:]))
        elif len(shape) == 4:  # first 2 dims are lists
            assert(len(ar) == shape[0]), \
                "Leading dim mismatch: %d != %d!" % (len(ar), shape[0])
            assert(shape[0] == 0 or len(ar[0]) == shape[1]), \
                "Second dim mismatch: %d != %d!" % (len(ar[0]), shape[1])
            assert(shape[0] == 0 or shape[1] == 0 or ar[0][0].shape == (shape[2], shape[3])), \
                "Shape mismatch: %s != %s!" % (str(ar[0][0].shape), str(shape[2:]))
        else:
            raise NotImplementedError("Number of dimensions must be <= 4!")


def lindblad_error_generators(dmbasis_ham, dmbasis_other, normalize,
                              other_mode="all"):
    """
    Compute the superoperator-generators corresponding to Lindblad terms.

    This routine computes the Hamiltonian and Non-Hamiltonian ("other")
    superoperator generators which correspond to the terms of the Lindblad
    expression:

    L(rho) = sum_i( h_i [A_i,rho] ) +
             sum_ij( o_ij * (B_i rho B_j^dag -
                             0.5( rho B_j^dag B_i + B_j^dag B_i rho) ) )

    where {A_i} and {B_i} are bases (possibly the same) for Hilbert Schmidt
    (density matrix) space with the identity element removed so that each
    A_i and B_i are traceless.  If we write L(rho) in terms of superoperators
    H_i and O_ij,

    L(rho) = sum_i( h_i H_i(rho) ) + sum_ij( o_ij O_ij(rho) )

    then this function computes the matrices for H_i and O_ij using the given
    density matrix basis.  Thus, if `dmbasis` is expressed in the standard
    basis (as it should be), the returned matrices are also in this basis.

    If these elements are used as projectors it may be usedful to normalize
    them (by setting `normalize=True`).  Note, however, that these projectors
    are not all orthogonal - in particular the O_ij's are not orthogonal to
    one another.

    Parameters
    ----------
    dmbasis_ham : list
        A list of basis matrices {B_i} *including* the identity as the first
        element, for the returned Hamiltonian-type error generators.  This
        argument is easily obtained by call to  :func:`pp_matrices` or a
        similar function.  The matrices are expected to be in the standard
        basis, and should be traceless except for the identity.  Matrices
        should be NumPy arrays or SciPy CSR sparse matrices.

    dmbasis_other : list
        A list of basis matrices {B_i} *including* the identity as the first
        element, for the returned Stochastic-type error generators.  This
        argument is easily obtained by call to  :func:`pp_matrices` or a
        similar function.  The matrices are expected to be in the standard
        basis, and should be traceless except for the identity.  Matrices
        should be NumPy arrays or SciPy CSR sparse matrices.

    normalize : bool
        Whether or not generators should be normalized so that
        numpy.linalg.norm(generator.flat) == 1.0  Note that the generators
        will still, in general, be non-orthogonal.

    other_mode : {"diagonal", "diag_affine", "all"}
        Which non-Hamiltonian Lindblad error generators to construct.
        Allowed values are: `"diagonal"` (only the diagonal Stochastic
        generators are returned; that is, the generators corresponding to the
        `i==j` terms in the Lindblad expression.), `"diag_affine"` (diagonal +
        affine generators), and `"all"` (all generators).


    Returns
    -------
    ham_generators : numpy.ndarray or list of SciPy CSR matrices
        If dense matrices where given, an array of shape (d-1,d,d), where d is
        the size of the basis, i.e. d == len(dmbasis).  `ham_generators[i]`
        gives the matrix for H_i.  If sparse matrices were given, a list
        of shape (d,d) CSR matrices.

    other_generators : numpy.ndarray or list of lists of SciPy CSR matrices
        If dense matrices where given, An array of shape (d-1,d-1,d,d),
        (2,d-1,d,d), or (d-1,d,d), where d is the size of the basis, for
        `other_mode` equal to `"all"`, `"diag_affine"`, or `"diagonal"`,
        respectively.  For instance, in the `"all"` case,
        `other_generators[i,j]` gives the matrix for O_ij.  If sparse matrices
        were given, the all but the final 2 dimensions are lists (e.g. the
        `"all"` case returns a list of lists of shape (d,d) CSR matrices).
    """
    if dmbasis_ham is not None:
        ham_mxs = dmbasis_ham  # list of basis matrices (assumed to be in std basis)
        ham_nMxs = len(ham_mxs)  # usually == d2, but not necessary (e.g. w/maxWeight)
    else:
        ham_nMxs = 0

    if dmbasis_other is not None:
        other_mxs = dmbasis_other  # list of basis matrices (assumed to be in std basis)
        other_nMxs = len(other_mxs)  # usually == d2, but not necessary (e.g. w/maxWeight)
    else:
        other_nMxs = 0

    if ham_nMxs > 0:
        d = ham_mxs[0].shape[0]
        sparse = _sps.issparse(ham_mxs[0])
    elif other_nMxs > 0:
        d = other_mxs[0].shape[0]
        sparse = _sps.issparse(other_mxs[0])
    else:
        d = 0  # will end up returning no generators
        sparse = False
    d2 = d**2
    normfn = _spsl.norm if sparse else _np.linalg.norm
    identityfn = (lambda d: _sps.identity(d, 'd', 'csr')) if sparse else _np.identity

    if ham_nMxs > 0 and other_nMxs > 0:
        assert(other_mxs[0].shape[0] == ham_mxs[0].shape[0]), \
            "Bases must have the same dimension!"

    if ham_nMxs > 0:
        assert(_np.isclose(normfn(ham_mxs[0] - identityfn(d) / _np.sqrt(d)), 0)),\
            "The first matrix in 'dmbasis_ham' must be the identity"

        hamLindbladTerms = [None] * (ham_nMxs - 1) if sparse else \
            _np.empty((ham_nMxs - 1, d2, d2), 'complex')

        for i, B in enumerate(ham_mxs[1:]):  # don't include identity
            hamLindbladTerms[i] = _lt.hamiltonian_to_lindbladian(B, sparse)  # in std basis
            if normalize:
                norm = normfn(hamLindbladTerms[i])  # same as norm(term.flat)
                if not _np.isclose(norm, 0):
                    hamLindbladTerms[i] /= norm  # normalize projector
                    assert(_np.isclose(normfn(hamLindbladTerms[i]), 1.0))
    else:
        hamLindbladTerms = None

    if other_nMxs > 0:
        assert(_np.isclose(normfn(other_mxs[0] - identityfn(d) / _np.sqrt(d)), 0)),\
            "The first matrix in 'dmbasis_other' must be the identity"

        if other_mode == "diagonal":
            otherLindbladTerms = [None] * (other_nMxs - 1) if sparse else \
                _np.empty((other_nMxs - 1, d2, d2), 'complex')
            for i, Lm in enumerate(other_mxs[1:]):  # don't include identity
                otherLindbladTerms[i] = _lt.nonham_lindbladian(Lm, Lm, sparse)
                if normalize:
                    norm = normfn(otherLindbladTerms[i])  # same as norm(term.flat)
                    if not _np.isclose(norm, 0):
                        otherLindbladTerms[i] /= norm  # normalize projector
                        assert(_np.isclose(normfn(otherLindbladTerms[i]), 1.0))

        elif other_mode == "diag_affine":
            otherLindbladTerms = [[None] * (other_nMxs - 1)] * 2 if sparse else \
                _np.empty((2, other_nMxs - 1, d2, d2), 'complex')
            for i, Lm in enumerate(other_mxs[1:]):  # don't include identity
                otherLindbladTerms[0][i] = _lt.nonham_lindbladian(Lm, Lm, sparse)
                otherLindbladTerms[1][i] = _lt.affine_lindbladian(Lm, sparse)
                if normalize:
                    for k in (0, 1):
                        norm = normfn(otherLindbladTerms[k][i])  # same as norm(term.flat)
                        if not _np.isclose(norm, 0):
                            otherLindbladTerms[k][i] /= norm  # normalize projector
                            assert(_np.isclose(normfn(otherLindbladTerms[k][i]), 1.0))

        else:  # other_mode == "all"
            otherLindbladTerms = \
                [[None] * (other_nMxs - 1) for i in range(other_nMxs - 1)] if sparse else \
                _np.empty((other_nMxs - 1, other_nMxs - 1, d2, d2), 'complex')

            for i, Lm in enumerate(other_mxs[1:]):  # don't include identity
                for j, Ln in enumerate(other_mxs[1:]):  # don't include identity
                    #print("DEBUG NONHAM LIND (%d,%d)" % (i,j)) #DEBUG!!!
                    otherLindbladTerms[i][j] = _lt.nonham_lindbladian(Lm, Ln, sparse)
                    if normalize:
                        norm = normfn(otherLindbladTerms[i][j])  # same as norm(term.flat)
                        if not _np.isclose(norm, 0):
                            otherLindbladTerms[i][j] /= norm  # normalize projector
                            assert(_np.isclose(normfn(otherLindbladTerms[i][j]), 1.0))
                    #I don't think this is true in general, but appears to be true for "pp" basis (why?)
                    #if j < i: # check that other[i,j] == other[j,i].C, i.e. other is Hermitian
                    #    assert(_np.isclose(_np.linalg.norm(
                    #                otherLindbladTerms[i][j]-
                    #                otherLindbladTerms[j][i].conjugate()),0))
    else:
        otherLindbladTerms = None

    #Check for orthogonality - otherLindblad terms are *not* orthogonal!
    #N = otherLindbladTerms.shape[0]
    #for i in range(N):
    #    for j in range(N):
    #        v1 = otherLindbladTerms[i,j].flatten()
    #        for k in range(N):
    #            for l in range(N):
    #                if k == i and l == j: continue
    #                v2 = otherLindbladTerms[k,l].flatten()
    #                if not _np.isclose(0, _np.vdot(v1,v2)):
    #                    print("%d,%d <-> %d,%d dot = %g [%g]" % (i,j,k,l,_np.vdot(v1,v2),_np.dot(v1,v2)))
    #                    #print("v1 = ",v1)
    #                    #print("v2 = ",v2)
    #                #    assert(False)
    #                #assert(_np.isclose(0, _np.vdot(v1,v2)))

    #Check hamiltonian error gens are orthogonal to others
    #N = otherLindbladTerms.shape[0]
    #for i,hlt in enumerate(hamLindbladTerms):
    #    v1 = hlt.flatten()
    #    for j in range(N):
    #        for k in range(N):
    #            v2 = otherLindbladTerms[j,k].flatten()
    #            assert(_np.isclose(0, _np.vdot(v1,v2)))

    return hamLindbladTerms, otherLindbladTerms


def lindblad_errgen_projections(errgen, ham_basis,
                                other_basis, mxBasis="gm",
                                normalize=True, return_generators=False,
                                other_mode="all", sparse=False):
    """
    Compute the projections of a gate error generator onto generators
    for the Lindblad-term errors when expressed in the given
    "projection basis".

    Parameters
    ----------
    errgen: : ndarray
      The error generator matrix to project.

    ham_basis: {'std', 'gm', 'pp', 'qt'}, list of matrices, or Basis object
        The basis 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.

    other_basis : {'std', 'gm', 'pp', 'qt'}, list of matrices, or Basis object
        The basis 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.

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

    normalize : bool, optional
      Whether or not the generators being projected onto are normalized, so
      that numpy.linalg.norm(generator.flat) == 1.0.  Note that the generators
      will still, in general, be non-orthogonal.

    return_generators : bool, optional
      If True, return the error generators projected against along with the
      projection values themseves.

    other_mode : {"diagonal", "diag_affine", "all"}
      Which non-Hamiltonian Lindblad error projections to obtain.
      Allowed values are: `"diagonal"` (only the diagonal Stochastic),
      `"diag_affine"` (diagonal + affine generators), and `"all"`
      (all generators).

    sparse : bool, optional
      Whether to create sparse or dense basis matrices when strings
      are given as `ham_basis` and `other_basis`

    Returns
    -------
    ham_projections : numpy.ndarray
      An array of length d-1, where d is the dimension of the gate,
      giving the projections onto the Hamiltonian-type Lindblad terms.

    other_projections : numpy.ndarray
      An array of shape (d-1,d-1), (2,d-1), or (d-1,), where d is the dimension
      of the gate, for `other_mode` equal to `"all"`, `"diag_affine"`, or
      `"diagonal"`, respectively.  Values give the projections onto the
      non-Hamiltonian-type Lindblad terms.

    ham_generators : numpy.ndarray
      The Hamiltonian-type Lindblad term generators, as would be returned
      from `lindblad_error_generators(pp_matrices(sqrt(d)), normalize)`.
      Shape is (d-1,d,d), and `ham_generators[i]` is in the standard basis.

    other_generators : numpy.ndarray
      The Stochastic-type Lindblad term generators, as would be returned
      from `lindblad_error_generators(pp_matrices(sqrt(d)), normalize)`.
      Shape is (d-1,d-1,d,d), (2,d-1,d,d), or (d-1,d,d) for `other_mode`
      equal to `"all"`, `"diag_affine"`, or `"diagonal"`, respectively,
      and `other_generators[i]` is in the std basis.

    """
    errgen_std = _bt.change_basis(errgen, mxBasis, "std")
    if _sps.issparse(errgen_std):
        errgen_std_flat = errgen_std.tolil().reshape(
            (errgen_std.shape[0] * errgen_std.shape[1], 1)).tocsr()  # b/c lil's are only type that can reshape...
    else:
        errgen_std_flat = errgen_std.flatten()
    errgen_std = None  # ununsed below, and sparse reshape doesn't copy, so mark as None

    d2 = errgen.shape[0]
    d = int(_np.sqrt(d2))
    #nQubits = _np.log2(d)

    #Get a list of the generators in corresspondence with the
    #  specified basis elements.
    if isinstance(ham_basis, _Basis):
        hamBasisMxs = ham_basis.elements
    elif isinstance(ham_basis, str):
        hamBasisMxs = _bt.basis_matrices(ham_basis, d2, sparse=sparse)
    else:
        hamBasisMxs = ham_basis

    if isinstance(other_basis, _Basis):
        otherBasisMxs = other_basis.elements
    elif isinstance(other_basis, str):
        otherBasisMxs = _bt.basis_matrices(other_basis, d2, sparse=sparse)
    else:
        otherBasisMxs = other_basis

    hamGens, otherGens = lindblad_error_generators(
        hamBasisMxs, otherBasisMxs, normalize, other_mode)  # in std basis

    if hamBasisMxs is not None:
        bsH = len(hamBasisMxs)  # basis size (not necessarily d2)
    else: bsH = 0

    if otherBasisMxs is not None:
        bsO = len(otherBasisMxs)  # basis size (not necessarily d2)
    else: bsO = 0

    if bsH > 0: sparse = _sps.issparse(hamBasisMxs[0])
    elif bsO > 0: sparse = _sps.issparse(otherBasisMxs[0])
    else: sparse = False  # default?

    assert(_np.isclose(d * d, d2))  # d2 must be a perfect square
    if bsH > 0:
        _assert_shape(hamGens, (bsH - 1, d2, d2), sparse)
    if bsO > 0:
        if other_mode == "diagonal":
            _assert_shape(otherGens, (bsO - 1, d2, d2), sparse)
        elif other_mode == "diag_affine":
            _assert_shape(otherGens, (2, bsO - 1, d2, d2), sparse)
        else:  # other_mode == "all"
            _assert_shape(otherGens, (bsO - 1, bsO - 1, d2, d2), sparse)

    #Perform linear least squares solve to find "projections" onto each otherGens element - defined so that
    #  sum_i projection_i * otherGen_i = (errgen_std-ham_errgen) as well as possible.

    #ham_error_gen = _np.einsum('i,ijk', hamProjs, hamGens)
    #other_errgen = errgen_std - ham_error_gen #what's left once hamiltonian errors are projected out

    #Do linear least squares soln to expressing errgen_std as a linear combo
    # of the lindblad generators
    if bsH > 0:
        if not sparse:
            H = hamGens.reshape((bsH - 1, d2**2)).T  # ham generators == columns
            Hdag = H.T.conjugate()

            #Do linear least squares: this is what takes the bulk of the time
            hamProjs = _np.linalg.solve(_np.dot(Hdag, H), _np.dot(Hdag, errgen_std_flat))
            hamProjs.shape = (hamGens.shape[0],)
        else:
            rows = [hamGen.tolil().reshape((1, d2**2)) for hamGen in hamGens]
            H = _sps.vstack(rows, 'csr').transpose()
            Hdag = H.copy().transpose().conjugate()

            #Do linear least squares: this is what takes the bulk of the time
            if _mt.safenorm(errgen_std_flat) < 1e-8:  # protect against singular RHS
                hamProjs = _np.zeros(bsH - 1, 'd')
            else:
                hamProjs = _spsl.spsolve(Hdag.dot(H), Hdag.dot(errgen_std_flat))
                if _sps.issparse(hamProjs): hamProjs = hamProjs.toarray().flatten()
            hamProjs.shape = (bsH - 1,)
    else:
        hamProjs = None

    if bsO > 0:
        if not sparse:
            if other_mode == "diagonal":
                O = otherGens.reshape((bsO - 1, d2**2)).T  # other generators == columns
            elif other_mode == "diag_affine":
                O = otherGens.reshape((2 * (bsO - 1), d2**2)).T  # other generators == columns
            else:
                O = otherGens.reshape(((bsO - 1)**2, d2**2)).T  # other generators == columns
            Odag = O.T.conjugate()

            #Do linear least squares: this is what takes the bulk of the time
            otherProjs = _np.linalg.solve(_np.dot(Odag, O), _np.dot(Odag, errgen_std_flat))

            if other_mode == "diagonal":
                otherProjs.shape = (otherGens.shape[0],)
            elif other_mode == "diag_affine":
                otherProjs.shape = (2, otherGens.shape[1])
            else:
                otherProjs.shape = (otherGens.shape[0], otherGens.shape[1])

        else:
            if other_mode == "diagonal":
                rows = [oGen.tolil().reshape((1, d2**2)) for oGen in otherGens]
                O = _sps.vstack(rows, 'csr').transpose()  # other generators == columns
            else:  # "diag_affine" or "all"
                rows = [oGen.tolil().reshape((1, d2**2)) for oGenRow in otherGens for oGen in oGenRow]
                O = _sps.vstack(rows, 'csr').transpose()  # other generators == columns
            Odag = O.copy().transpose().conjugate()  # TODO: maybe conjugate copies data?

            #Do linear least squares: this is what takes the bulk of the time
            if _mt.safenorm(errgen_std_flat) < 1e-8:  # protect against singular RHS
                if other_mode == "diagonal": otherProjs = _np.zeros(bsO - 1, 'd')
                elif other_mode == "diag_affine": otherProjs = _np.zeros((2, bsO - 1), 'd')
                else: otherProjs = _np.zeros((bsO - 1, bsO - 1), 'd')
            else:
                otherProjs = _spsl.spsolve(Odag.dot(O), Odag.dot(errgen_std_flat))
                if _sps.issparse(otherProjs): otherProjs = otherProjs.toarray().flatten()

            if other_mode == "diagonal":
                otherProjs.shape = (bsO - 1,)
            elif other_mode == "diag_affine":
                otherProjs.shape = (2, bsO - 1)
            else:  # other_mode == "all"
                otherProjs.shape = (bsO - 1, bsO - 1)
    else:
        otherProjs = None

    #check err gens are linearly independent -- but can take a very long time, so comment out!
    #assert(_np.linalg.matrix_rank(H,1e-7) == H.shape[1])
    #assert(_np.linalg.matrix_rank(O,1e-7) == O.shape[1])
    #if False: # further check against older (slower) version
    #    M = _np.concatenate( (hamGens.reshape((bs-1,d2**2)).T, otherGens.reshape(((bs-1)**2,d2**2)).T), axis=1)
    #    assert(_np.linalg.matrix_rank(M,1e-7) == M.shape[1]) #check err gens are linearly independent
    #    Mdag = M.T.conjugate()
    #    print("DB D: %.1f" % (time.time()-t)); t = time.time()
    #    projs = _np.linalg.solve(_np.dot(Mdag,M), _np.dot(Mdag,errgen_std_flat))
    #    hamProjs_chk = projs[0:(bs-1)]
    #    otherProjs_chk = projs[(bs-1):]
    #    assert(_np.linalg.norm(hamProjs-hamProjs_chk) < 1e-6)
    #    assert(_np.linalg.norm(otherProjs-otherProjs_chk) < 1e-6)

    if return_generators:
        return hamProjs, otherProjs, hamGens, otherGens
    else:
        return hamProjs, otherProjs


def projections_to_lindblad_terms(hamProjs, otherProjs, ham_basis, other_basis,
                                  other_mode="all", return_basis=True):
    """
    Converts the projections of an error generator onto basis elements into
    the Lindblad-term dictionary and basis used to individually specify
    Lindblad terms.

    Parameters
    ----------
    hamProjs : numpy.ndarray
        An array of length d-1, where d is the dimension of the projected error
        generator, giving the projections onto the Hamiltonian-type Lindblad
        terms.

    otherProjs : numpy.ndarray
        An array of shape (d-1,d-1), (2,d-1), or (d-1,), where d is the dimension
        of the projected error generator, for `other_mode` equal to `"all"`,
        `"diag_affine"`, or `"diagonal"`, respectively.  Values give the
        projections onto the non-Hamiltonian-type Lindblad terms.

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

    other_basis : {'std', 'gm', 'pp', 'qt'}, list of matrices, or Basis object
        The basis used to construct `otherProjs`.  Allowed values are
        Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt),
        list of numpy arrays, or a custom basis object.

    other_mode : {"diagonal", "diag_affine", "all"}
        Which non-Hamiltonian Lindblad error projections `otherProjs` includes.
        Allowed values are: `"diagonal"` (only the diagonal Stochastic),
        `"diag_affine"` (diagonal + affine generators), and `"all"`
        (all generators).

    return_basis : bool, optional
        Whether to return a :class:`Basis` containing the elements
        corresponding to labels within the returned `Ltermdict`.


    Returns
    -------
    Ltermdict : dict
        Keys are `(termType, basisLabel1, <basisLabel2>)`
        tuples, where `termType` is `"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 have 1 basis label
        to indicate a *diagonal* term and otherwise have 2 basis labels to
        specify off-diagonal non-Hamiltonian Lindblad terms.  Basis labels
        are taken from `ham_basis` and `other_basis`.  Values are complex
        coefficients (the projections).

    basis : Basis
        A single basis containing all the basis labels used in `Ltermdict` (and
        *only* those elements).  Only returned when `return_basis == True`.
    """
    assert(not (ham_basis is None and other_basis is None)), \
        "At least one of `ham_basis` and `other_basis` must be non-None"

    # Make None => length-0 arrays so iteration code works below (when basis is None)
    if hamProjs is None: hamProjs = _np.empty(0, 'd')
    if otherProjs is None:
        otherProjs = _np.empty(0, 'd') if other_mode == "diagonal" \
            else _np.empty((0, 0), 'd')

    # Construct a pair of dictionaries describing all of the
    # Lindblad-terms:
    #   Ltermdict keys= ('H',basisLbl), ('S',basisLbl), or ('S',bLbl1,bLbl2)
    #             vals= coefficients of these terms (projections from errgen)
    #   basisdict keys= basis labels (just has to match Ltermdict keys)
    #             vals= basis matrices - can be either sparse or dense
    Ltermdict = _collections.OrderedDict()
    basisdict = _collections.OrderedDict()

    if return_basis:
        def set_basis_el(blbl, bel):
            """ Sets an elment of basisdict, checking for consistency """
            if blbl in basisdict:
                assert(_mt.safenorm(basisdict[blbl] - bel) < 1e-8), "Ambiguous basis el label %s" % blbl
            else:
                basisdict[blbl] = bel
    else:
        def set_basis_el(blbl, bel):
            pass

    #Add Hamiltonian error elements
    if ham_basis is not None:
        ham_lbls = ham_basis.labels
        ham_mxs = ham_basis.elements  # can be sparse
        assert(len(ham_mxs[1:]) == len(hamProjs))
        for coeff, lbl, bmx in zip(hamProjs, ham_lbls[1:], ham_mxs[1:]):  # skip identity
            Ltermdict[('H', lbl)] = coeff
            set_basis_el(lbl, bmx)
    else:
        ham_lbls = []

    #Add "other" error elements
    if other_basis is not None:
        other_lbls = other_basis.labels
        other_mxs = other_basis.elements  # can be sparse
        if other_mode == "diagonal":
            assert(len(other_mxs[1:]) == len(otherProjs))
            for coeff, lbl, bmx in zip(otherProjs, other_lbls[1:], other_mxs[1:]):  # skip identity
                Ltermdict[('S', lbl)] = coeff
                set_basis_el(lbl, bmx)

        elif other_mode == "diag_affine":
            assert((2, len(other_mxs[1:])) == otherProjs.shape)
            for coeff, lbl, bmx in zip(otherProjs[0], other_lbls[1:], other_mxs[1:]):  # skip identity
                Ltermdict[('S', lbl)] = coeff
                set_basis_el(lbl, bmx)
            for coeff, lbl, bmx in zip(otherProjs[1], other_lbls[1:], other_mxs[1:]):  # skip identity
                Ltermdict[('A', lbl)] = coeff
                set_basis_el(lbl, bmx)

        else:
            assert((len(other_mxs[1:]), len(other_mxs[1:])) == otherProjs.shape)
            for i, (lbl1, bmx1) in enumerate(zip(other_lbls[1:], other_mxs[1:])):  # skip identity
                set_basis_el(lbl1, bmx1)
                for j, (lbl2, bmx2) in enumerate(zip(other_lbls[1:], other_mxs[1:])):  # skip identity
                    set_basis_el(lbl2, bmx2)
                    Ltermdict[('S', lbl1, lbl2)] = otherProjs[i, j]
    else:
        other_lbls = []

    #Turn basisdict into a Basis to return
    if return_basis:
        if ham_basis == other_basis:
            basis = ham_basis
        elif ham_basis is None or set(ham_lbls).issubset(set(other_lbls)):
            basis = other_basis
        elif other_basis is None or set(other_lbls).issubset(set(ham_lbls)):
            basis = ham_basis
        else:
            #Create an ExplictBasis using the matrices in basisdict plus the identity
            sparse = True; real = True
            if ham_basis is not None:
                elshape = ham_basis.elshape
                sparse = sparse and ham_basis.sparse
                real = real and ham_basis.real
            if other_basis is not None:
                elshape = other_basis.elshape
                sparse = sparse and other_basis.sparse
                real = real and other_basis.real

            d = elshape[0]
            Id = _sps.identity(d, 'complex', 'csr') / _np.sqrt(d) if sparse \
                else _np.identity(d, 'complex') / _np.sqrt(d)

            lbls = ['I'] + list(basisdict.keys())
            mxs = [Id] + list(basisdict.values())
            basis = _ExplicitBasis(mxs, lbls, name=None,
                                   real=real, sparse=sparse)
        return Ltermdict, basis
    else:
        return Ltermdict


def lindblad_terms_to_projections(Ltermdict, basis, other_mode="all"):
    """
    Convert a set of Lindblad terms into a dense matrix/grid of projections.

    Essentially the inverse of :function:`projections_to_lindblad_terms`.

    Parameters
    ----------
    Ltermdict : dict
        A dictionary specifying which Linblad terms are present in the gate
        parameteriztion.  Keys are `(termType, basisLabel1, <basisLabel2>)`
        tuples, where `termType` is `"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).

    basis : Basis, optional
        A basis mapping the labels used in the keys of `Ltermdict` to
        basis matrices (e.g. numpy arrays or Scipy sparse matrices).  The
        first element of this basis should be an identity element, and
        will be propagated to the returned `ham_basis` and `other_basis`.

    other_mode : {"diagonal", "diag_affine", "all"}
      Which non-Hamiltonian terms are allowed in `Ltermdict`.
      Allowed values are: `"diagonal"` (only the diagonal Stochastic),
      `"diag_affine"` (diagonal + affine generators), and `"all"`
      (all generators).

    Returns
    -------
    hamProjs : numpy.ndarray
        An array of length `basisdim-1`, giving the projections onto a
        full set of the Hamiltonian-type Lindblad terms (onto each element of
        `ham_basis`).

    otherProjs : numpy.ndarray
        An array of shape (d-1,d-1), (2,d-1), or (d-1,), where d=`basisdim`
        for `other_mode` equal to `"all"`, `"diag_affine"`, or `"diagonal"`,
        respectively.  Values give the projections onto the non-Hamiltonian
        -type Lindblad terms.

    ham_basis: Basis
        The basis used to construct `hamProjs`.

    other_basis : Basis
        The basis used to construct `otherProjs`.

    hamBasisIndices : OrderedDict
        A dictionary mapping the some or all of the basis labels of `basisdict`
        to the integers 0 to `len(ham_basis)`.  These are indices into
        `hamProjs`, giving the projection associated with each Hamiltonian
        basis element.

    otherBasisIndices : OrderedDict
        A dictionary mapping the some or all of the basis labels of `basisdict`
        to the integers 0 to `len(other_basis)`.  These are row and column
        indices into `otherProjs`, giving the projection associated with each
        pair of "other" basis elements (or single basis element if
        `other_mode!="all"`).
    """
    #Separately enumerate the (distinct) basis elements used for Hamiltonian
    # and non-Hamiltonian error terms
    #print("DB: lindblad term to proj: \n",Ltermdict,"\n",basis)
    hamBasisLabels = []
    otherBasisLabels = []
    for termLbl, coeff in Ltermdict.items():
        if isinstance(termLbl, str): termLbl = (termLbl[0], termLbl[1:])  # e.g. "HXX" => ('H','XX')
        termType = termLbl[0]
        if termType == "H":  # Hamiltonian
            assert(len(termLbl) == 2), "Hamiltonian term labels should have form ('H',<basis element label>)"
            if termLbl[1] not in hamBasisLabels:
                hamBasisLabels.append(termLbl[1])

        elif termType == "S":  # Stochastic
            if other_mode in ("diagonal", "diag_affine"):
                assert(len(termLbl) == 2), "Stochastic term labels should have form ('S',<basis element label>)"
                if termLbl[1] not in otherBasisLabels:
                    otherBasisLabels.append(termLbl[1])
            else:
                assert(len(termLbl) == 3), "Stochastic term labels should have form ('S',<bel1>, <bel2>)"
                if termLbl[1] not in otherBasisLabels:
                    otherBasisLabels.append(termLbl[1])
                if termLbl[2] not in otherBasisLabels:
                    otherBasisLabels.append(termLbl[2])

        elif termType == "A":  # Affine
            assert(other_mode == "diag_affine"), "Affine labels are only allowed in an affine mode"
            assert(len(termLbl) == 2), "Affine term labels should have form ('A',<basis element label>)"
            if termLbl[1] not in otherBasisLabels:
                otherBasisLabels.append(termLbl[1])

    #Construct bases
    # Note: the lists of basis matrices shouldn't contain the identity, since
    # the terms above shouldn't contain identity terms - but `basis` should
    # contain an identity element as it's first element, so add this identity el
    # to non-empty bases (empty bases stay empty!) to be consistent with the
    # rest of the framework (bases *have* Ids)

    sparse = basis.sparse
    if set(hamBasisLabels) == set(basis.labels):
        ham_basis = basis
    else:
        Id = basis[0]
        ham_basis_mxs = [basis[bl] for bl in hamBasisLabels]
        if len(ham_basis_mxs) > 0:
            ham_basis = _ExplicitBasis([Id] + ham_basis_mxs, ['I'] + hamBasisLabels,
                                       name=None, real=True, sparse=sparse)
        else:
            ham_basis = _ExplicitBasis(ham_basis_mxs, name=None, real=True, sparse=sparse)

    if set(otherBasisLabels) == set(basis.labels):
        other_basis = basis
    else:
        Id = basis[0]
        other_basis_mxs = [basis[bl] for bl in otherBasisLabels]
        if len(other_basis_mxs) > 0:
            other_basis = _ExplicitBasis([Id] + other_basis_mxs, ['I'] + otherBasisLabels,
                                         name=None, real=True, sparse=sparse)
        else:
            other_basis = _ExplicitBasis(other_basis_mxs, name=None, real=True, sparse=sparse)

    bsH, bsO = len(ham_basis), len(other_basis)
    #print("DB: constructed ham_basis = ",ham_basis)
    #print("DB: other basis = ",other_basis)

    #Create projection (term coefficient) arrays - or return None if
    # the corresponding basis is empty (as per our convention)
    hamProjs = _np.zeros(bsH - 1, 'complex') if bsH > 0 else None
    if bsO > 0:
        if other_mode == "diagonal":  # OK if this runs for 'auto' too since then len(otherBasisIndices) == 0
            otherProjs = _np.zeros(bsO - 1, 'complex')
        elif other_mode == "diag_affine":
            otherProjs = _np.zeros((2, bsO - 1), 'complex')
        else:
            otherProjs = _np.zeros((bsO - 1, bsO - 1), 'complex')
    else: otherProjs = None

    #Fill arrays
    hamBasisIndices = {lbl: i - 1 for i, lbl in enumerate(ham_basis.labels)}      # -1 to compensate for identity as
    otherBasisIndices = {lbl: i - 1 for i, lbl in enumerate(other_basis.labels)}  # first element (not in projections).
    for termLbl, coeff in Ltermdict.items():
        if isinstance(termLbl, str): termLbl = (termLbl[0], termLbl[1:])  # e.g. "HXX" => ('H','XX')
        termType = termLbl[0]
        if termType == "H":  # Hamiltonian
            k = hamBasisIndices[termLbl[1]]  # index of coefficient in array
            hamProjs[k] = coeff
        elif termType == "S":  # Stochastic
            if other_mode == "diagonal":
                k = otherBasisIndices[termLbl[1]]  # index of coefficient in array
                otherProjs[k] = coeff
            elif other_mode == "diag_affine":
                k = otherBasisIndices[termLbl[1]]  # index of coefficient in array
                otherProjs[0, k] = coeff
            else:  # other_mode == "all"
                k = otherBasisIndices[termLbl[1]]  # index of row in "other" coefficient matrix
                j = otherBasisIndices[termLbl[2]]  # index of col in "other" coefficient matrix
                otherProjs[k, j] = coeff
        elif termType == "A":  # Affine
            assert(other_mode == "diag_affine")
            k = otherBasisIndices[termLbl[1]]  # index of coefficient in array
            otherProjs[1, k] = coeff

    return hamProjs, otherProjs, ham_basis, other_basis


def lindblad_projections_to_paramvals(hamProjs, otherProjs, param_mode="cptp",
                                      other_mode="all", truncate=True):
    """
    Construct the array of Lindblad-gate parameter values from the separate
    arrays of Hamiltonian and non-Hamiltonian Lindblad-term projections.

    When `cptp=True`, this function handles parameterizing the projections
    to that for (real) parameter values correspond to projections for a valid
    CPTP gate (e.g. by parameterizing the Cholesky decomposition of `otherProjs`
    instead of otherProjs itself).  This function is closely related to
    implementation details of the LindbladOp class.

    Parameters
    ----------
    hamProjs : numpy.ndarray
        An array of length d-1, where d is the gate dimension, giving the
        projections onto a full set of the Hamiltonian-type Lindblad terms.

    otherProjs : numpy.ndarray
        An array of shape (d-1,d-1), (2,d-1), or (d-1,), where d is the gate
        dimension, for `other_mode` equal to `"all"`,`"diag_affine"`, or
        `"diagonal"`, respectively.  Values give the projections onto a full
        set of non-Hamiltonian-type Lindblad terms.

    param_mode : {"unconstrained", "cptp", "depol", "reldepol"}
        Describes how values in `hamProjs` and `otherProj` relate to the
        returned parameter values.  Allowed values are:
        `"unconstrained"` (projs are independent unconstrained parameters),
        `"cptp"` (independent parameters but constrained so map is CPTP),
        `"reldepol"` (all non-Ham. diagonal projs take the *same* value),
        `"depol"` (same as `"reldepol"` but projs must be *positive*)

    other_mode : {"diagonal", "diag_affine", "all"}
        Which non-Hamiltonian Lindblad error projections `otherProjs` includes.
        Allowed values are: `"diagonal"` (only the diagonal Stochastic),
        `"diag_affine"` (diagonal + affine generators), and `"all"`.

    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 projections
        cannot be parameterized as specified.

    Returns
    -------
    numpy.ndarray
        A 1D array of real parameter values consisting of d-1 Hamiltonian
        values followed by either (d-1)^2, 2*(d-1), or just d-1 non-Hamiltonian
        values for `other_mode` equal to `"all"`, `"diag_affine"`, or
        `"diagonal"`, respectively.
    """
    if hamProjs is not None:
        assert(_np.isclose(_np.linalg.norm(hamProjs.imag), 0)), \
            "Hamiltoian projections (coefficients) are not all real!"
        hamParams = hamProjs.real
    else:
        hamParams = _np.empty(0, 'd')

    if otherProjs is not None:
        if other_mode == "diagonal":
            assert(_np.isclose(_np.linalg.norm(_np.imag(otherProjs)), 0)), \
                "Diagonal stochastic projections (coefficients) are not all real!"

            if param_mode == "depol":  # otherParams is a *single-element* 1D vector of the sqrt of each diagonal el
                assert(truncate or all([v >= -1e-12 for v in otherProjs])), \
                    "Lindblad coefficients are not CPTP (truncate == False)!"
                assert(truncate or all([_np.isclose(v, otherProjs[0]) for v in otherProjs])), \
                    "Diagonal lindblad coefficients are not equal (truncate == False)!"
                otherProj = _np.mean(otherProjs.clip(1e-16, 1e100))
                otherParams = _np.array(_np.sqrt(_np.real(otherProj)), 'd')  # shape (1,)

            elif param_mode == "cptp":  # otherParams is a 1D vector of the sqrts of diagonal els
                assert(truncate or all([v >= -1e-12 for v in otherProjs])), \
                    "Lindblad coefficients are not CPTP (truncate == False)!"
                otherProjs = otherProjs.clip(1e-16, 1e100)
                otherParams = _np.sqrt(otherProjs.real)  # shape (bsO-1,)
            else:  # "unconstrained": otherParams is a 1D vector of the real diagonal els of otherProjs
                otherParams = otherProjs.real  # shape (bsO-1,)

        elif other_mode == "diag_affine":
            assert(_np.isclose(_np.linalg.norm(_np.imag(otherProjs)), 0)), \
                "Diagonal stochastic and affine projections (coefficients) are not all real!"

            if param_mode == "depol":  # otherParams is a single depol value + unconstrained affine coeffs
                assert(truncate or all([v >= -1e-12 for v in otherProjs[0]])), \
                    "Lindblad coefficients are not CPTP (truncate == False)!"
                assert(truncate or all([_np.isclose(v, otherProjs[0, 0]) for v in otherProjs[0]])), \
                    "Diagonal lindblad coefficients are not equal (truncate == False)!"
                depolProj = _np.mean(otherProjs[0, :].clip(1e-16, 1e100))
                otherParams = _np.concatenate(([_np.sqrt(_np.real(depolProj))],
                                               otherProjs[1].real))  # shape (1+(bsO-1),)

            elif param_mode == "cptp":  # Note: does not constrained affine coeffs to CPTP
                assert(truncate or all([v >= -1e-12 for v in otherProjs[0]])), \
                    "Lindblad coefficients are not CPTP (truncate == False)!"
                diagParams = _np.sqrt(_np.real(otherProjs[0, :]).clip(1e-16, 1e100))  # shape (bsO-1,)
                otherParams = _np.concatenate((diagParams, otherProjs[1].real))  # diag + affine params

            else:  # param_mode == "unconstrained": otherParams is a 1D vector of the real diagonal els of otherProjs
                otherParams = otherProjs.real  # shape (2,bsO-1)

        else:  # other_mode == "all"
            assert(_np.isclose(_np.linalg.norm(otherProjs - otherProjs.T.conjugate()), 0)
                   ), "Other projection/coefficient mx is not Hermitian!"
            assert(param_mode != "depol"), "`depol` is not supported when `other_mode == 'all'`"

            bsO = otherProjs.shape[0] + 1  # +1 to keep convention that this is the basis (w/Identity) size
            otherParams = _np.empty((bsO - 1, bsO - 1), 'd')

            if param_mode == "cptp":  # otherParams mx stores Cholesky decomp

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

                assert(truncate or all([ev >= -1e-12 for ev in evals])), \
                    "Lindblad coefficients are not CPTP (truncate == False)!"

                pos_evals = evals.clip(1e-16, 1e100)
                otherProjs = _np.dot(U, _np.dot(_np.diag(pos_evals), Ui))
                try:
                    Lmx = _np.linalg.cholesky(otherProjs)

                # if Lmx not postitive definite, try again with 1e-12 (same lines as above)
                except _np.linalg.LinAlgError:                         # pragma: no cover
                    pos_evals = evals.clip(1e-12, 1e100)                # pragma: no cover
                    otherProjs = _np.dot(U, _np.dot(_np.diag(pos_evals), Ui))  # pragma: no cover
                    Lmx = _np.linalg.cholesky(otherProjs)                  # pragma: no cover

                for i in range(bsO - 1):
                    assert(_np.linalg.norm(_np.imag(Lmx[i, i])) < IMAG_TOL)
                    otherParams[i, i] = Lmx[i, i].real
                    for j in range(i):
                        otherParams[i, j] = Lmx[i, j].real
                        otherParams[j, i] = Lmx[i, j].imag

            else:  # param_mode == "unconstrained": otherParams mx stores otherProjs (hermitian) directly
                for i in range(bsO - 1):
                    assert(_np.linalg.norm(_np.imag(otherProjs[i, i])) < IMAG_TOL)
                    otherParams[i, i] = otherProjs[i, i].real
                    for j in range(i):
                        otherParams[i, j] = otherProjs[i, j].real
                        otherParams[j, i] = otherProjs[i, j].imag
    else:
        otherParams = _np.empty(0, 'd')

    assert(not _np.iscomplexobj(hamParams))   # params should always
    assert(not _np.iscomplexobj(otherParams))  # be *real*
    return _np.concatenate((hamParams, otherParams.flat))


def paramvals_to_lindblad_projections(paramvals, ham_basis_size,
                                      other_basis_size, param_mode="cptp",
                                      other_mode="all", Lmx=None):
    """
    Construct the separate arrays of Hamiltonian and non-Hamiltonian
    Lindblad-term projections from the array of Lindblad-gate parameter values.

    This function essentially performs the inverse of
    :function:`lindblad_projections_to_paramvals`.

    Parameters
    ----------
    paramvals : numpy.ndarray
        A 1D array of real parameter values consisting of d-1 Hamiltonian
        values followed by either (d-1)^2 or just d-1 non-Hamiltonian
        values (the latter when `other_mode in ('diagonal','diag_affine')`).

    ham_basis_size, other_basis_size : int
        The number of elements in the Hamiltonian and non-Hamiltonian
        bases used to construct `paramvals`.  As such, `ham_basis_size`
        gives the offset into `paramvals` where the non-Hamiltonian
        parameters begin.

    param_mode : {"unconstrained", "cptp", "depol", "reldepol"}
        Specifies how the Lindblad-term coefficients are mapped to the set of
        (real) parameter values.  This really just applies to the "other"
        (non-Hamiltonian) coefficients.  "unconstrained" means that ranging
        over the parameter values lets the coefficient-matrix vary over all
        matrices, "cptp" restricts this to postitive matrices. "depol"
        maps all of the coefficients to the *same, positive* parameter (only
        available for "diagonal" and "diag_affine" other-modes), and "reldepol"
        does the same thing but without the positivity constraint.

    other_mode : {"all", "diagonal", "diag_affine"}
        Specifies the structure of the matrix of other (non-Hamiltonian)
        coefficients.  If d is the gate dimension, "all" means a (d-1,d-1)
        matrix is used; "diagonal" means just the (d2-1,) diagonal of this
        matrix is used; "diag_affine" means the coefficients are in a (2,d2-1)
        array with the diagonal-term coefficients being the first row and the
        affine coefficients being the second row.

    Lmx : ndarray, optional
        Scratch space that is used to store the lower-triangular
        Cholesky decomposition matrix that is used to construct
        the "other" projections when there is a CPTP constraint.

    Returns
    -------
    hamProjs : numpy.ndarray
        An array of length d-1, where d is the gate dimension, giving the
        projections onto a full set of the Hamiltonian-type Lindblad terms.

    otherProjs : numpy.ndarray
        An array of shape (d-1,d-1) or (d-1,) or (2,d-1) where d is the gate
        dimension, giving the projections onto a full set of non-Hamiltonian
        -type Lindblad terms (see `other_mode` above).
    """
    bsH = ham_basis_size
    bsO = other_basis_size

    if Lmx is None:
        Lmx = _np.zeros((bsO - 1, bsO - 1), 'complex') if bsO > 0 else None

    # self.paramvals = [hamCoeffs] + [otherParams]
    #  where hamCoeffs are *real* and of length d2-1 (self.dim == d2)
    if bsH > 0:
        hamCoeffs = paramvals[0:bsH - 1]
        nHam = bsH - 1
    else:
        hamCoeffs = None
        nHam = 0

    #built up otherCoeffs based on param_mode and nonham_mode
    if bsO > 0:
        if other_mode == "diagonal":
            otherParams = paramvals[nHam:]
            expected_shape = (1,) if (param_mode in ("depol", "reldepol")) else (bsO - 1,)
            assert(otherParams.shape == expected_shape)
            if param_mode in ("depol", "reldepol"):
                otherParams = otherParams[0] * _np.ones(bsO - 1, 'd')  # replicate single param bsO-1 times

            if param_mode in ("cptp", "depol"):
                otherCoeffs = otherParams**2  # Analagous to L*L_dagger
            else:  # "unconstrained"
                otherCoeffs = otherParams

        elif other_mode == "diag_affine":

            if param_mode in ("depol", "reldepol"):
                otherParams = paramvals[nHam:].reshape((1 + bsO - 1,))
                otherCoeffs = _np.empty((2, bsO - 1), 'd')  # leave as real type b/c doesn't have complex entries
                if param_mode == "depol":
                    otherCoeffs[0, :] = otherParams[0]**2
                else:
                    otherCoeffs[0, :] = otherParams[0]
                otherCoeffs[1, :] = otherParams[1:]

            else:
                otherParams = paramvals[nHam:].reshape((2, bsO - 1))
                if param_mode == "cptp":
                    otherCoeffs = otherParams.copy()
                    otherCoeffs[0, :] = otherParams[0]**2
                else:  # param_mode == "unconstrained"
                    #otherCoeffs = _np.empty((2,bsO-1),'complex')
                    otherCoeffs = otherParams

        else:  # other_mode == "all"
            otherParams = paramvals[nHam:].reshape((bsO - 1, bsO - 1))

            if param_mode == "cptp":
                #  otherParams is an array of length (bs-1)*(bs-1) that
                #  encodes a lower-triangular matrix "Lmx" via:
                #  Lmx[i,i] = otherParams[i,i]
                #  Lmx[i,j] = otherParams[i,j] + 1j*otherParams[j,i] (i > j)
                for i in range(bsO - 1):
                    Lmx[i, i] = otherParams[i, i]
                    for j in range(i):
                        Lmx[i, j] = otherParams[i, j] + 1j * otherParams[j, i]

                #The matrix of (complex) "other"-coefficients is build by
                # assuming Lmx is its Cholesky decomp; means otherCoeffs
                # is pos-def.

                # NOTE that the Cholesky decomp with all positive real diagonal
                # elements is *unique* for a given positive-definite otherCoeffs
                # matrix, but we don't care about this uniqueness criteria and so
                # the diagonal els of Lmx can be negative and that's fine -
                # otherCoeffs will still be posdef.
                otherCoeffs = _np.dot(Lmx, Lmx.T.conjugate())

                #DEBUG - test for pos-def
                #evals = _np.linalg.eigvalsh(otherCoeffs)
                #DEBUG_TOL = 1e-16; #print("EVALS DEBUG = ",evals)
                #assert(all([ev >= -DEBUG_TOL for ev in evals]))

            else:  # param_mode == "unconstrained"
                #otherParams holds otherCoeff real and imaginary parts directly
                otherCoeffs = _np.empty((bsO - 1, bsO - 1), 'complex')
                for i in range(bsO - 1):
                    otherCoeffs[i, i] = otherParams[i, i]
                    for j in range(i):
                        otherCoeffs[i, j] = otherParams[i, j] + 1j * otherParams[j, i]
                        otherCoeffs[j, i] = otherParams[i, j] - 1j * otherParams[j, i]
    else:
        otherCoeffs = None

    return hamCoeffs, otherCoeffs


#TODO: replace two_qubit_gate, one_qubit_gate, unitary_to_pauligate_* with
# calls to this one and unitary_to_processmx
def rotation_gate_mx(r, mxBasis="gm"):
    """
    Construct a rotation operation matrix.

    Build the operation matrix corresponding to the unitary
    `exp(-i * (r[0]/2*PP[0]*sqrt(d) + r[1]/2*PP[1]*sqrt(d) + ...) )`
    where `PP' is the array of Pauli-product matrices
    obtained via `pp_matrices(d)`, where `d = sqrt(len(r)+1)`.
    The division by 2 is for convention, and the sqrt(d) is to
    essentially un-normalise the matrices returned by `pp_matrices`
    to they are equal to products of the *standard* Pauli matrices.

    Parameters
    ----------
    r : tuple
        A tuple of coeffiecients, one per non-identity
        Pauli-product basis element

    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).
.
    Returns
    -------
    numpy array
        a d^2 x d^2 operation matrix in the specified basis.
    """
    d = int(round(_np.sqrt(len(r) + 1)))
    assert(d**2 == len(r) + 1), "Invalid number of rotation angles"

    #get Pauli-product matrices (in std basis)
    pp = _bt.basis_matrices('pp', d**2)
    assert(len(r) == len(pp[1:]))

    #build unitary (in std basis)
    ex = _np.zeros((d, d), 'complex')
    for rot, pp_mx in zip(r, pp[1:]):
        ex += rot / 2.0 * pp_mx * _np.sqrt(d)
    U = _spl.expm(-1j * ex)
    stdGate = unitary_to_process_mx(U)

    ret = _bt.change_basis(stdGate, 'std', mxBasis)

    return ret


def project_model(model, targetModel,
                  projectiontypes=('H', 'S', 'H+S', 'LND'),
                  genType="logG-logT"):
    """
    Construct one or more new models by projecting the error generator of
    `model` onto some sub-space then reconstructing.

    Parameters
    ----------
    model : Model
        The model whose error generator should be projected.

    targetModel : Model
        The set of target (ideal) gates.

    projectiontypes : tuple of {'H','S','H+S','LND','LNDCP'}
        Which projections to use.  The length of this tuple gives the
        number of `Model` objects returned.  Allowed values are:

        - 'H' = Hamiltonian errors
        - 'S' = Stochastic Pauli-channel errors
        - 'H+S' = both of the above error types
        - 'LND' = errgen projected to a normal (CPTP) Lindbladian
        - 'LNDF' = errgen projected to an unrestricted (full) Lindbladian

    genType : {"logG-logT", "logTiG"}
      The type of error generator to compute.  Allowed values are:

      - "logG-logT" : errgen = log(gate) - log(target_op)
      - "logTiG" : errgen = log( dot(inv(target_op), gate) )

    Returns
    -------
    projected_models : list of Models
       Elements are projected versions of `model` corresponding to
       the elements of `projectiontypes`.

    Nps : list of parameter counts
       Integer parameter counts for each model in `projected_models`.
       Useful for computing the expected log-likelihood or chi2.
    """

    opLabels = list(model.operations.keys())  # operation labels
    basis = model.basis

    #The projection basis needs to be a basis for density matrices
    # (i.e. 2x2 mxs in 1Q case) rather than superoperators (4x4 mxs
    # in 1Q case) - whcih is what model.basis is.  So, we just extract
    # a builtin basis name for the projection basis.
    if basis.name in ('pp', 'gm', 'std', 'qt'):
        proj_basis_name = basis.name
    else:
        proj_basis_name = 'pp'  # model.basis is weird so just use paulis as projection basis

    if basis.name != targetModel.basis.name:
        raise ValueError("Basis mismatch between model (%s) and target (%s)!"
                         % (model.basis.name, targetModel.basis.name))

    # Note: set to "full" parameterization so we can set the gates below
    #  regardless of what parameterization the original model had.
    gsDict = {}; NpDict = {}
    for p in projectiontypes:
        gsDict[p] = model.copy()
        gsDict[p].set_all_parameterizations("full")
        NpDict[p] = 0

    errgens = [error_generator(model.operations[gl],
                               targetModel.operations[gl],
                               targetModel.basis, genType)
               for gl in opLabels]

    for gl, errgen in zip(opLabels, errgens):
        if ('H' in projectiontypes) or ('H+S' in projectiontypes):
            hamProj, hamGens = std_errgen_projections(
                errgen, "hamiltonian", proj_basis_name, basis, True)
            #ham_error_gen = _np.einsum('i,ijk', hamProj, hamGens)
            ham_error_gen = _np.tensordot(hamProj, hamGens, (0, 0))
            ham_error_gen = _bt.change_basis(ham_error_gen, "std", basis)

        if ('S' in projectiontypes) or ('H+S' in projectiontypes):
            stoProj, stoGens = std_errgen_projections(
                errgen, "stochastic", proj_basis_name, basis, True)
            #sto_error_gen = _np.einsum('i,ijk', stoProj, stoGens)
            sto_error_gen = _np.tensordot(stoProj, stoGens, (0, 0))
            sto_error_gen = _bt.change_basis(sto_error_gen, "std", basis)

        if ('LND' in projectiontypes) or ('LNDF' in projectiontypes):
            HProj, OProj, HGens, OGens = \
                lindblad_errgen_projections(
                    errgen, proj_basis_name, proj_basis_name, basis, normalize=False,
                    return_generators=True)
            #Note: return values *can* be None if an empty/None basis is given
            #lnd_error_gen = _np.einsum('i,ijk', HProj, HGens) + \
            #                _np.einsum('ij,ijkl', OProj, OGens)
            lnd_error_gen = _np.tensordot(HProj, HGens, (0, 0)) + \
                _np.tensordot(OProj, OGens, ((0, 1), (0, 1)))
            lnd_error_gen = _bt.change_basis(lnd_error_gen, "std", basis)

        targetOp = targetModel.operations[gl]

        if 'H' in projectiontypes:
            gsDict['H'].operations[gl] = operation_from_error_generator(
                ham_error_gen, targetOp, genType)
            NpDict['H'] += len(hamProj)

        if 'S' in projectiontypes:
            gsDict['S'].operations[gl] = operation_from_error_generator(
                sto_error_gen, targetOp, genType)
            NpDict['S'] += len(stoProj)

        if 'H+S' in projectiontypes:
            gsDict['H+S'].operations[gl] = operation_from_error_generator(
                ham_error_gen + sto_error_gen, targetOp, genType)
            NpDict['H+S'] += len(hamProj) + len(stoProj)

        if 'LNDF' in projectiontypes:
            gsDict['LNDF'].operations[gl] = operation_from_error_generator(
                lnd_error_gen, targetOp, genType)
            NpDict['LNDF'] += HProj.size + OProj.size

        if 'LND' in projectiontypes:
            evals, U = _np.linalg.eig(OProj)
            pos_evals = evals.clip(0, 1e100)  # clip negative eigenvalues to 0
            OProj_cp = _np.dot(U, _np.dot(_np.diag(pos_evals), _np.linalg.inv(U)))
            #OProj_cp is now a pos-def matrix
            #lnd_error_gen_cp = _np.einsum('i,ijk', HProj, HGens) + \
            #                   _np.einsum('ij,ijkl', OProj_cp, OGens)
            lnd_error_gen_cp = _np.tensordot(HProj, HGens, (0, 0)) + \
                _np.tensordot(OProj_cp, OGens, ((0, 1), (0, 1)))
            lnd_error_gen_cp = _bt.change_basis(lnd_error_gen_cp, "std", basis)

            gsDict['LND'].operations[gl] = operation_from_error_generator(
                lnd_error_gen_cp, targetOp, genType)
            NpDict['LND'] += HProj.size + OProj.size

        #Removed attempt to contract H+S to CPTP by removing positive stochastic projections,
        # but this doesn't always return the gate to being CPTP (maybe b/c of normalization)...
        #sto_error_gen_cp = _np.einsum('i,ijk', stoProj.clip(None,0), stoGens)
        #  # (only negative stochastic projections OK)
        #sto_error_gen_cp = _tools.std_to_pp(sto_error_gen_cp)
        #gsHSCP.operations[gl] = _tools.operation_from_error_generator(
        #    ham_error_gen, targetOp, genType) #+sto_error_gen_cp

    #DEBUG!!!
    #print("DEBUG: BEST sum neg evals = ",_tools.sum_of_negative_choi_evals(model))
    #print("DEBUG: LNDCP sum neg evals = ",_tools.sum_of_negative_choi_evals(gsDict['LND']))

    #Check for CPTP where expected
    #assert(_tools.sum_of_negative_choi_evals(gsHSCP) < 1e-6)
    #assert(_tools.sum_of_negative_choi_evals(gsDict['LND']) < 1e-6)

    #Collect and return requrested results:
    ret_gs = [gsDict[p] for p in projectiontypes]
    ret_Nps = [NpDict[p] for p in projectiontypes]
    return ret_gs, ret_Nps


def get_a_best_case_gauge_transform(gate_mx, target_gate_mx, returnAll=False):
    """
    Returns a gauge transformation that maps `gate_mx` into a matrix that is
    co-diagonal with `target_gate_mx`, i.e. they share a common set of eigenvectors.

    Gauge transformations effectively change the basis of all the gates in a model.
    From the perspective of a single gate a gauge transformation leaves it's
    eigenvalues the same and changes its eigenvectors.  This function finds a *real*
    transformation that transforms the eigenspaces of `gate_mx` so that there exists
    a set of eigenvectors which diagonalize both `gate_mx` and `target_gate_mx`.

    Parameters
    ----------
    gate_mx, target_gate_mx : numpy.ndarray
        The gate and target-gate matrices.

    returnAll : bool, optional
        If true, also return the matrices of eigenvectors
        for `Ugate` for gate_mx and `Utgt` for target_gate_mx such
        that `U = dot(Utgt, inv(Ugate))` is real.

    Returns
    -------
    U : numpy.ndarray
        A gauge transformation such that if `epgate = U * gate_mx * U_inv`,
        then `epgate` (which has the same eigenalues as `gate_mx`), can be
        diagonalized with a set of eigenvectors that also diagonalize
        `target_gate_mx`.  Furthermore, `U` is real.

    Ugate, Utgt : numpy.ndarray
        only if `returnAll == True`.  See above.
    """

    # A complication that must be dealt with is that
    # the eigenvalues of `target_gate_mx` can be degenerate,
    # and so matching up eigenvalues can't be done *just* based on value.
    # Our algorithm consists of two steps:
    # 1) match gate & target eigenvalues based on value, ensuring conjugacy
    #    relationships between eigenvalues are preserved.
    # 2) for each eigenvalue/vector of `gate`, project the eigenvector onto
    #    the eigenspace of `tgt_gate` corresponding to the matched eigenvalue.
    #    (treat conj-pair eigenvalues of `gate` together).

    # we want a matrix that gauge-transforms gate_mx into a matrix as
    # close to target_gate_mx as possible, i.e. that puts gate_mx's
    # eigenvalues in the eigenspaces of target_gate_mx.  This is done
    # by Ubest = _np.dot(Utgt, inv(Uop)), but there are often degrees
    # of freedom in Uop because of its degeneracies.  Also, we want Ubest
    # to be *real*, so we need to ensure the conjugacy structure of Utgt
    # and Uop match...

    assert(_np.linalg.norm(gate_mx.imag) < 1e-8)
    assert(_np.linalg.norm(target_gate_mx.imag) < 1e-8)

    if True:  # NEW approach that gives sorted eigenvectors
        def get_eigenspace_pairs(mx, TOL=1e-6):
            evals, U = _np.linalg.eig(mx)  # so mx = U * evals * Uinv
            espace_pairs = {}; conj_pair_indices = []

            #Pass 1: real evals and positive-imaginary-element-of-conjugate pair evals
            #  (these are the representatives of "eigenspace pairs")
            for i, ev in enumerate(evals):
                if ev.imag < -TOL:
                    conj_pair_indices.append(i); continue  # save for pass2

                #see if ev is already in espace_pairs
                for k, v in espace_pairs.items():
                    if abs(k - ev) < TOL:
                        espace_pairs[k]['indices'].append(i)
                        espace_pairs[k]['conj_pair_indices'].append(None)
                        #espace_pairs[k]['evecs'].append(U[:,i])
                        break
                else:
                    espace_pairs[ev] = {'indices': [i], 'conj_pair_indices': [None]}

            #Pass 2: negative-imaginary-part elements of evals that occur in conjugate pairs
            for i in conj_pair_indices:
                ev_pos = _np.conjugate(evals[i])
                for k, v in espace_pairs.items():  # ev_pos *should* be in espace_pairs
                    if abs(k - ev_pos) < TOL:
                        #found the correct eigenspace-pair to add this eval & evec to,
                        # now figure our where to put this index based on conjugacy relationships,
                        # i.e. U[:,esp['indices'][i]] is always conjugate to U[:,esp['conj_pair_indices'][i]]
                        for jj, j in enumerate(espace_pairs[k]['indices']):
                            if espace_pairs[k]['conj_pair_indices'][jj] is None:  # an empty slot
                                espace_pairs[k]['conj_pair_indices'][jj] = i
                                U[:, i] = U[:, j].conj()
                                break
                        else:
                            raise ValueError("Nowhere to place a conjugate eigenvector %d-dim eigenbasis for %s!"
                                             % (len(espace_pairs[k]['indices']), str(k)))

                        break
                else:
                    raise ValueError("Expected to find %s as an espace-pair representative in %s"
                                     % (str(ev_pos), str(espace_pairs.keys())))

            #if not (_np.allclose(mx, _np.dot(U, _np.dot(_np.diag(evals), _np.linalg.inv(U))))):
            #    import bpdb; bpdb.set_trace()
            return evals, U, espace_pairs

        def standard_diag(mx, TOL=1e-6):
            evals, U, espairs = get_eigenspace_pairs(mx)
            std_evals = []
            std_evecs = []
            sorted_rep_evals = sorted(list(espairs.keys()), key=lambda x: (x.real, x.imag))
            for ev in sorted_rep_evals:  # iterate in sorted order just for definitiveness
                info = espairs[ev]
                dim = len(info['indices'])  # dimension of this eigenspace (and it's pair, if there is one)

                #Ensure real eigenvalue blocks should have real eigenvectors
                if abs(ev.imag) < TOL:
                    #find linear combinations of the eigenvectors that are real
                    Usub = U[:, info['indices']]
                    if _np.linalg.norm(Usub.imag) > TOL:
                        # Im part of Usub * combo = Usub.real*combo.imag + Usub.imag*combo.real
                        combo_real_imag = _mt.nullspace(_np.concatenate((Usub.imag, Usub.real), axis=1))
                        combos = combo_real_imag[0:dim, :] + 1j * combo_real_imag[dim:, :]
                        if combos.shape[1] != dim:
                            raise ValueError(("Can only find %d (< %d) *real* linear combinations of"
                                              " vectors in eigenspace for %s!") % (combos.shape[1], dim, str(ev)))
                        U[:, info['indices']] = _np.dot(Usub, combos)
                        assert(_np.linalg.norm(U[:, info['indices']].imag) < TOL)

                    #Add real eigenvalues and vectors
                    std_evals.extend([ev] * dim)
                    std_evecs.extend([U[:, i] for i in info['indices']])

                else:  # complex eigenvalue case - should have conjugate pair info
                    #Ensure blocks for conjugate-pairs of eigenvalues follow one after another and
                    # corresponding eigenvectors (e.g. the first of each block) are conjugate pairs
                    # (this is already done in the eigenspace construction)
                    assert(len(info['conj_pair_indices']) == dim)
                    std_evals.extend([ev] * dim)
                    std_evals.extend([_np.conjugate(ev)] * dim)
                    std_evecs.extend([U[:, i] for i in info['indices']])
                    std_evecs.extend([U[:, i] for i in info['conj_pair_indices']])

            return _np.array(std_evals), _np.array(std_evecs).T

        #Create "gate_tilde" which has the eigenvectors of gate_mx around the matched eigenvalues of target_gate_mx
        # Doing this essentially decouples the problem of eigenvalue matching from the rest of the task -
        # after gate_tilde is created, it and target_gate_mx have exactly the *same* eigenvalues.
        evals_tgt, Utgt = _np.linalg.eig(target_gate_mx)
        evals_gate, Uop = _np.linalg.eig(gate_mx)
        pairs = _mt.minweight_match_realmxeigs(evals_gate, evals_tgt)
        replace_evals = _np.array([evals_tgt[j] for _, j in pairs])
        gate_tilde = _np.dot(Uop, _np.dot(_np.diag(replace_evals), _np.linalg.inv(Uop)))

        #Create "standard diagonalizations" of gate_tilde and target_gate_mx, which give
        # sort the eigenvalues and ensure eigenvectors occur in *corresponding* conjugate pairs
        # (e.g. even when evals +1j and -1j have multiplicity 4, the first 4-D eigenspace, the
        evals_tgt, Utgt = standard_diag(target_gate_mx)
        evals_tilde, Uop = standard_diag(gate_tilde)
        assert(_np.allclose(evals_tgt, evals_tilde))

        #Update Utgt so that Utgt * inv_Uop is close to the identity
        kite = _mt.get_kite(evals_tgt)  # evals are grouped by standard_diag, so this works
        D_prior_to_proj = _np.dot(_np.linalg.inv(Utgt), Uop)
        #print("D prior to projection to ",kite," kite:"); _mt.print_mx(D_prior_to_proj)
        D = _mt.project_onto_kite(D_prior_to_proj, kite)
        start = 0
        for i, k in enumerate(kite):
            slc = slice(start, start + k)
            dstart = start + k
            for kk in kite[i + 1:]:
                if k == kk and _np.isclose(evals_tgt[start], evals_tgt[dstart].conj()):  # conjugate block!
                    dslc = slice(dstart, dstart + kk)
                    # enforce block conjugacy needed to retain Uproj conjugacy structure
                    D[dslc, dslc] = D[slc, slc].conj()
                    break
                dstart += kk
            start += k
        Utgt = _np.dot(Utgt, D)  # update Utgt

        Utrans = _np.dot(Utgt, _np.linalg.inv(Uop))
        assert(_np.linalg.norm(_np.imag(Utrans)) < 1e-7)
        Utrans = Utrans.real  # _np.real_if_close(Utrans, tol=1000)

        if returnAll:
            return Utrans, Uop, Utgt, evals_tgt
        else:
            return Utrans

    evals_tgt, Utgt = _np.linalg.eig(target_gate_mx)
    evals_gate, Uop = _np.linalg.eig(gate_mx)

    #_, pairs = _mt.minweight_match(evals_tgt, evals_gate, return_pairs=True)
    pairs = _mt.minweight_match_realmxeigs(evals_tgt, evals_gate)

    #Form eigenspaces of Utgt
    eigenspace = {}  # key = index of target eigenval, val = assoc. eigenspace
    for i, ev in enumerate(evals_tgt):
        for j in eigenspace:
            if _np.isclose(ev, evals_tgt[j]):  # then add evector[i] to this eigenspace
                eigenspace[j].append(Utgt[:, i])
                eigenspace[i] = eigenspace[j]  # reference!
                break
        else:
            eigenspace[i] = [Utgt[:, i]]  # new list = new eigenspace

    #Project each eigenvector (col of Uop) onto space of cols
    evectors = {}  # key = index of gate eigenval, val = assoc. (projected) eigenvec
    for ipair, (i, j) in enumerate(pairs):
        #print("processing pair (i,j) = ",i,j)
        if j in evectors: continue  # we already processed this one!

        # non-orthog projection:
        # v = E * coeffs s.t. |E*coeffs-v|^2 is minimal  (E is not square so can't invert)
        # --> E.dag * v = E.dag * E * coeffs
        # --> inv(E.dag * E) * E.dag * v = coeffs
        # E*coeffs = E * inv(E.dag * E) * E.dag * v
        E = _np.array(eigenspace[i]).T; Edag = E.T.conjugate()
        coeffs = _np.dot(_np.dot(_np.linalg.inv(_np.dot(Edag, E)), Edag), Uop[:, j])
        evectors[j] = _np.dot(E, coeffs)

        #check for conjugate pair
        #DB: print("Looking for conjugate:")
        for i2, j2 in pairs[ipair + 1:]:
            if abs(evals_gate[j].imag) > 1e-6 and _np.isclose(evals_gate[j], _np.conjugate(evals_gate[j2])) \
               and _np.allclose(Uop[:, j], Uop[:, j2].conj()):
                #DB: print("Found conjugate at j = ",j2)
                evectors[j2] = _np.conjugate(evectors[j])
                # x = _np.linalg.solve(_np.dot(Edag, E), _np.dot(Edag, evectors[j2]))
                #assert(_np.isclose(_np.linalg.norm(x),_np.linalg.norm(coeffs))) ??
                #check that this vector is in the span of eigenspace[i2]?

    #build new "Utgt" using specially chosen linear combos of degenerate-eigenvecs
    Uproj = _np.array([evectors[i] for i in range(Utgt.shape[1])]).T
    assert(_np.allclose(_np.dot(Uproj, _np.dot(_np.diag(evals_tgt), _np.linalg.inv(Uproj))), target_gate_mx))

    #This is how you get the eigenspace-projected gate:
    #  epgate = _np.dot(Uproj, _np.dot(_np.diag(evals_gate), Uproj_inv))
    #  epgate = _np.real_if_close(epgate, tol=1000)

    # G = Uop * evals_gate * Uop_inv  => eval_gate = Uop_inv * G * Uop
    # epgate = Uproj * evals_gate * Uproj_inv  (eigenspace-projected gate)
    # so  epgate = (Uproj Uop_inv) G (Uproj Uop_inv)_inv => (Uproj Uop_inv) is
    # a "best_gauge_transform" for G, i.e. it makes G codiagonal with G_tgt
    Ubest = _np.dot(Uproj, _np.linalg.inv(Uop))
    assert(_np.linalg.norm(_np.imag(Ubest)) < 1e-7)
    # this should never happen & indicates an uncaught failure in
    # minweight_match_realmxeigs(...)

    Ubest = Ubest.real

    if returnAll:
        return Ubest, Uop, Uproj, evals_tgt
    else:
        return Ubest


def project_to_target_eigenspace(model, targetModel, EPS=1e-6):
    """
    Project each gate of `model` onto the eigenspace of the corresponding
    gate within `targetModel`.  Return the resulting `Model`.

    Parameters
    ----------
    model, targetModel : Model
        The model being projected and the model specifying the "target"
        eigen-spaces, respectively.

    EPS : float, optional
        Small magnitude specifying how much to "nudge" the target gates
        before eigen-decomposing them, so that their spectra will have the
        same conjugacy structure as the gates of `model`.

    Returns
    -------
    Model
    """
    ret = targetModel.copy()
    ret.set_all_parameterizations("full")  # so we can freely assign gates new values

    for gl, gate in model.operations.items():
        tgt_gate = targetModel.operations[gl]

        #Essentially, we want to replace the eigenvalues of `tgt_gate`
        # (and *only* the eigenvalues) with those of `gate`.  This is what
        # a "best gate gauge transform does" (by definition)
        gate_mx = gate.todense()
        Ugauge = get_a_best_case_gauge_transform(gate_mx, tgt_gate.todense())
        Ugauge_inv = _np.linalg.inv(Ugauge)

        epgate = _np.dot(Ugauge, _np.dot(gate_mx, Ugauge_inv))
        ret.operations[gl] = epgate

    return ret


def unitary_to_pauligate(U):
    """
    Get the linear operator on (vectorized) density
    matrices corresponding to a n-qubit unitary
    operator on states.

    Parameters
    ----------
    U : numpy array
        A dxd array giving the action of the unitary
        on a state in the sigma-z basis.
        where d = 2 ** n-qubits

    Returns
    -------
    numpy array
        The operator on density matrices that have been
        vectorized as d**2 vectors in the Pauli basis.
    """
    assert U.shape[0] == U.shape[1], '"Unitary" matrix is not square'
    return _bt.change_basis(unitary_to_process_mx(U), 'std', 'pp')


def is_valid_lindblad_paramtype(typ):
    """
    Whether `typ` is a recognized Lindblad-gate parameterization type.

    A *Lindblad type* is comprised of a parameter specification followed
    optionally by an evolution-type suffix.  The parameter spec can be
    "GLND" (general unconstrained Lindbladian), "CPTP" (cptp-constrained),
    or any/all of the letters "H" (Hamiltonian), "S" (Stochastic, CPTP),
    "s" (Stochastic), "A" (Affine), "D" (Depolarization, CPTP),
    "d" (Depolarization) joined with plus (+) signs.  Note that "H"
    cannot appear alone, and that "A" cannot appear without one of
    {"S","s","D","d"}. The suffix can be non-existent (density-matrix),
    "terms" (state-vector terms) or "clifford terms" (stabilizer-state
    terms).  For example, valid Lindblad types are "H+S", "H+d+A",
    "CPTP clifford terms", or "S+A terms".

    Returns
    -------
    bool
    """
    try:
        baseTyp, _ = split_lindblad_paramtype(typ)
    except ValueError:
        return False  # if can't even split `typ`
    return baseTyp in ("CPTP", "H+S", "S", "H+S+A", "S+A", "H+D", "D", "H+D+A", "D+A",
                       "GLND", "H+s", "s", "H+s+A", "s+A", "H+d", "d", "H+d+A", "d+A")


def split_lindblad_paramtype(typ):
    """
    Splits a Lindblad-gate parameteriation type into
    a base-type (e.g. "H+S") and an evolution-type
    string.

    Parameters
    ----------
    typ : str
        The parameterization type, e.g. "H+S terms".

    Returns
    -------
    base_type : str
        The "base-parameterization" part of `typ`.
    evotype : str
        The evolution type corresponding to `typ`.
    """
    bTyp = typ.split()[0]  # "base" type
    evostr = " ".join(typ.split()[1:])

    if evostr == "": evotype = "densitymx"
    elif evostr == "terms": evotype = "svterm"
    elif evostr == "clifford terms": evotype = "cterm"
    else: raise ValueError("Unrecognized evotype in `typ`=%s" % typ)
    return bTyp, evotype


def eLabelToOutcome(povm_and_effect_lbl):
    """TODO: Docstring """
    # Helper fn: POVM_ELbl:sslbls -> Elbl mapping
    if povm_and_effect_lbl is None:
        return "NONE"  # Dummy label for placeholding
    else:
        if isinstance(povm_and_effect_lbl, _Label):
            last_underscore = povm_and_effect_lbl.name.rindex('_')
            effect_lbl = povm_and_effect_lbl.name[last_underscore + 1:]
        else:
            last_underscore = povm_and_effect_lbl.rindex('_')
            effect_lbl = povm_and_effect_lbl[last_underscore + 1:]
        return effect_lbl  # effect label alone *is* the outcome


def eLabelToPOVM(povm_and_effect_lbl):
    """TODO: Docstring """
    # Helper fn: POVM_ELbl:sslbls -> POVM mapping
    if povm_and_effect_lbl is None:
        return "NONE"  # Dummy label for placeholding
    else:
        if isinstance(povm_and_effect_lbl, _Label):
            last_underscore = povm_and_effect_lbl.name.rindex('_')
            povm_name = povm_and_effect_lbl.name[:last_underscore]
        else:
            last_underscore = povm_and_effect_lbl.rindex('_')
            povm_name = povm_and_effect_lbl[:last_underscore]
        return povm_name