def _steadystate_svd_dense(L, ss_args):
    """
    Find the steady state(s) of an open quantum system by solving for the
    nullspace of the Liouvillian.
    """
    atol = 1e-12
    rtol = 1e-12
    if settings.debug:
        print('Starting SVD solver...')

    u, s, vh = svd(L.full(), full_matrices=False)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T

    if ss_args['all_states']:
        rhoss_list = []
        for n in range(ns.shape[1]):
            rhoss = Qobj(vec2mat(ns[:, n]), dims=L.dims[0])
            rhoss_list.append(rhoss / rhoss.tr())
        return rhoss_list

    else:
        rhoss = Qobj(vec2mat(ns[:, 0]), dims=L.dims[0])
        return rhoss / rhoss.tr()

def propagator_steadystate(U):
    """Find the steady state for successive applications of the propagator
    :math:`U`.

    Parameters
    ----------
    U : qobj
        Operator representing the propagator.

    Returns
    -------
    a : qobj
        Instance representing the steady-state density matrix.

    """

    evals, evecs = la.eig(U.full())

    ev_min, ev_idx = _get_min_and_index(abs(evals - 1.0))

    evecs = evecs.T
    rho = Qobj(vec2mat(evecs[ev_idx]), dims=U.dims[0])
    rho = rho * (1.0 / rho.tr())
    rho = 0.5 * (rho + rho.dag())  # make sure rho is herm
    return rho

def _steadystate_eigen(L, ss_args):
    """
    Internal function for solving the steady state problem by
    finding the eigenvector corresponding to the zero eigenvalue
    of the Liouvillian using ARPACK.
    """
    if settings.debug:
        print('Starting Eigen solver...')

    dims = L.dims[0]
    shape = prod(dims[0])
    L = L.data.tocsc()

    if ss_args['use_rcm']:
        if settings.debug:
            old_band = sp_bandwidth(L)[0]
            print('Original bandwidth:', old_band)
        perm = reverse_cuthill_mckee(L)
        rev_perm = np.argsort(perm)
        L = sp_permute(L, perm, perm, 'csc')
        if settings.debug:
            rcm_band = sp_bandwidth(L)[0]
            print('RCM bandwidth:', rcm_band)
            print('Bandwidth reduction factor:', round(old_band/rcm_band, 1))

    eigval, eigvec = eigs(L, k=1, sigma=1e-15, tol=ss_args['tol'],
                          which='LM', maxiter=ss_args['maxiter'])

    if ss_args['use_rcm']:
        eigvec = eigvec[np.ix_(rev_perm,)]

    data = vec2mat(eigvec)
    data = 0.5 * (data + data.conj().T)
    out = Qobj(data, dims=dims, isherm=True)
    return out/out.tr()

def test_QobjDiagonals():
    "Qobj diagonals"
    data = np.random.random(
        (5, 5)) + 1j * np.random.random((5, 5)) - (0.5 + 0.5j)
    A = Qobj(data)
    b = A.diag()
    assert_(np.all(b - np.diag(data) == 0))

def _init_ptrace_stuff(sel):
    psi0 = Qobj(config.psi0,
                dims=config.psi0_dims,
                shape=config.psi0_shape)
    qtf90.qutraj_run.init_ptrace_stuff(config.psi0_dims[0],
                                       np.array(sel) + 1,
                                       psi0.ptrace(sel).shape[0])

def floquet_master_equation_tensor(Alist, f_energies):
    """
    Construct a tensor that represents the master equation in the floquet
    basis (with constant Hamiltonian and collapse operators).

    Simplest RWA approximation [Grifoni et al, Phys.Rep. 304 229 (1998)]

    Parameters
    ----------

    Alist : list
        A list of Floquet-Markov master equation rate matrices.

    f_energies : array
        The Floquet energies.

    Returns
    -------

    output : array

        The Floquet-Markov master equation tensor `R`.

    """

    if isinstance(Alist, list):
        # Alist can be a list of rate matrices corresponding
        # to different operators that couple to the environment
        N, M = np.shape(Alist[0])
    else:
        # or a simple rate matrix, in which case we put it in a list
        Alist = [Alist]
        N, M = np.shape(Alist[0])

    R = Qobj(scipy.sparse.csr_matrix((N * N, N * N)), [[N, N], [N, N]],
             [N * N, N * N])

    R.data = R.data.tolil()
    for I in range(N * N):
        a, b = vec2mat_index(N, I)
        for J in range(N * N):
            c, d = vec2mat_index(N, J)

            R.data[I, J] = -1.0j * (f_energies[a] - f_energies[b]) * \
                (a == c) * (b == d)

            for A in Alist:
                s1 = s2 = 0
                for n in range(N):
                    s1 += A[a, n] * (n == c) * (n == d) - A[n, a] * \
                        (a == c) * (a == d)
                    s2 += (A[n, a] + A[n, b]) * (a == c) * (b == d)

                dR = (a == b) * s1 - 0.5 * (1 - (a == b)) * s2

                if dR != 0.0:
                    R.data[I, J] += dR

    R.data = R.data.tocsr()
    return R

def test_QobjExpmExplicitlySparse():
    "Qobj expm (sparse)"
    data = np.random.random(
        (15, 15)) + 1j * np.random.random((15, 15)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.expm(method='sparse')
    assert_((B.data.todense() - np.matrix(la.expm(data)) < 1e-10).all())

def rand_ket(N, density=1, dims=None):
    """Creates a random Nx1 sparse ket vector.

    Parameters
    ----------
    N : int
        Number of rows for output quantum operator.
    density : float
        Density between [0,1] of output ket state.
    dims : list
        Left-dimensions of quantum object.  Used for specifying
        tensor structure. Default is dims=[[N]].

    Returns
    -------
    oper : qobj
        Nx1 ket state quantum operator.

    """
    if dims:
        _check_ket_dims(dims, N)
    X = sp.rand(N, 1, density, format="csr")
    X.data = X.data - 0.5
    Y = X.copy()
    Y.data = 1.0j * (np.random.random(len(X.data)) - 0.5)
    X = X + Y
    X.sort_indices()
    X = Qobj(X)
    if dims:
        return Qobj(X / X.norm(), dims=[dims, [1]], shape=[N, 1])
    else:
        return Qobj(X / X.norm())

def _steadystate_svd_dense(L, ss_args):
    """
    Find the steady state(s) of an open quantum system by solving for the
    nullspace of the Liouvillian.
    """
    ss_args['info'].pop('weight', None)
    atol = 1e-12
    rtol = 1e-12
    if settings.debug:
        logger.debug('Starting SVD solver.')
    _svd_start = time.time()
    u, s, vh = svd(L.full(), full_matrices=False)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    _svd_end = time.time()
    ss_args['info']['solution_time'] = _svd_end-_svd_start
    if ss_args['all_states']:
        rhoss_list = []
        for n in range(ns.shape[1]):
            rhoss = Qobj(vec2mat(ns[:, n]), dims=L.dims[0])
            rhoss_list.append(rhoss / rhoss.tr())
        if ss_args['return_info']:
            return rhoss_list, ss_args['info']
        else:
            if ss_args['return_info']:
                return rhoss_list, ss_args['info']
            else:
                return rhoss_list
    else:
        rhoss = Qobj(vec2mat(ns[:, 0]), dims=L.dims[0])
        return rhoss / rhoss.tr()

def test_QobjMulNonsquareDims():
    """
    Qobj: multiplication w/ non-square qobj.dims

    Checks for regression of #331.
    """
    data = np.array([[0, 1], [1, 0]])

    q1 = Qobj(data)
    q1.dims[0].append(1)
    q2 = Qobj(data)

    assert_equal((q1 * q2).dims, [[2, 1], [2]])
    assert_equal((q2 * q1.dag()).dims, [[2], [2, 1]])

    # Note that this is [[2], [2]] instead of [[2, 1], [2, 1]],
    # as matching dimensions of 1 are implicitly partial traced out.
    # (See #331.)
    assert_equal((q1 * q2 * q1.dag()).dims, [[2], [2]])
    
    # Because of the above, we also need to check for extra indices
    # that aren't of length 1.
    q1 = Qobj([[ 1.+0.j,  0.+0.j],
         [ 0.+0.j,  1.+0.j],
         [ 0.+0.j,  1.+0.j],
         [ 1.+0.j,  0.+0.j],
         [ 0.+0.j,  0.-1.j],
         [ 0.+1.j,  0.+0.j],
         [ 1.+0.j,  0.+0.j],
         [ 0.+0.j, -1.+0.j]
     ], dims=[[4, 2], [2]])

    assert_equal((q1 * q2 * q1.dag()).dims, [[4, 2], [4, 2]])

def test_Qobj_sqrtm():
    "Qobj sqrtm"
    data = np.random.random(
        (5, 5)) + 1j * np.random.random((5, 5)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.sqrtm()
    assert_(A == B * B)

def test_QobjHerm():
    "Qobj Hermicity"
    N = 10
    data = np.random.random(
        (N, N)) + 1j * np.random.random((N, N)) - (0.5 + 0.5j)
    q = Qobj(data)
    assert_equal(q.isherm, False)

    data = data + data.conj().T
    q = Qobj(data)
    assert_(q.isherm)

    q_a = destroy(5)
    assert_(not q_a.isherm)

    q_ad = create(5)
    assert_(not q_ad.isherm)

    # test addition of two nonhermitian operators adding up to a hermitian one
    q_x = q_a + q_ad
    assert_(q_x.isherm)  # isherm use the _isherm cache from q_a + q_ad
    q_x._isherm = None   # reset _isherm cache
    assert_(q_x.isherm)  # recalculate _isherm

    # test addition of one hermitan and one nonhermitian operator
    q = q_x + q_a
    assert_(not q.isherm)
    q._isherm = None
    assert_(not q.isherm)

    # test addition of two hermitan operators
    q = q_x + q_x
    assert_(q.isherm)
    q._isherm = None
    assert_(q.isherm)

def test_QobjFull():
    "Qobj full"
    data = np.random.random(
        (15, 15)) + 1j * np.random.random((15, 15)) - (0.5 + 0.5j)
    A = Qobj(data)
    b = A.full()
    assert_(np.all(b - data == 0))

def test_QobjExpm():
    "Qobj expm (dense)"
    data = np.random.random(
        (15, 15)) + 1j * np.random.random((15, 15)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.expm()
    assert_((B.data.todense() - np.matrix(la.expm(data)) < 1e-10).all())

def _steadystate_svd_dense(L, atol=1e-12, rtol=0, all_steadystates=False,
                           verbose=False):
    """
    Find the steady state(s) of an open quantum system by solving for the
    nullspace of the Liouvillian.
    """
    if verbose:
        print('Starting SVD solver...')
        start_time = time.time()

    u, s, vh = svd(L.full(), full_matrices=False)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T

    if verbose:
        print('SVD solver time: ', time.time() - start_time)

    if all_steadystates:
        rhoss_list = []
        for n in range(ns.shape[1]):
            rhoss = Qobj(vec2mat(ns[:, n]), dims=L.dims[0])
            rhoss_list.append(rhoss / rhoss.tr())
        return rhoss_list

    else:
        rhoss = Qobj(vec2mat(ns[:, 0]), dims=L.dims[0])
        return rhoss / rhoss.tr()

def operator_to_vector(op):
    """
    Create a vector representation of a quantum operator given
    the matrix representation.
    """
    q = Qobj()
    q.dims = [op.dims, [1]]
    q.data = sp_reshape(op.data.T, (np.prod(op.shape), 1))
    return q

def vector_to_operator(op):
    """
    Create a matrix representation given a quantum operator in
    vector form.
    """
    q = Qobj()
    q.dims = op.dims[0]
    q.data = sp_reshape(op.data.T, q.shape).T
    return q

def test_QobjDagger():
    "Qobj adjoint (dagger)"
    data = np.random.random((5, 5)) + 1j * np.random.random((5, 5)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.dag()
    assert_(np.all(B.data.todense() - np.matrix(data.conj().T) == 0))
    assert_equal(A.isherm, B.isherm)
    assert_equal(A.type, B.type)
    assert_equal(A.superrep, B.superrep)

def test_QobjConjugate():
    "Qobj conjugate"
    data = np.random.random((5, 5)) + 1j * np.random.random((5, 5)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.conj()
    assert_(np.all(B.data.todense() - np.matrix(data.conj()) == 0))
    assert_equal(A.isherm, B.isherm)
    assert_equal(A.type, B.type)
    assert_equal(A.superrep, B.superrep)

def vector_to_operator(op):
    """
    Create a matrix representation given a quantum operator in
    vector form.
    """
    q = Qobj()
    q.dims = op.dims[0]
    n = int(np.sqrt(op.shape[0]))
    q.data = sp_reshape(op.data.T, (n, n)).T
    return q