r"""Return a Hamilton operator ``H`` and a minimal list of collapse
    operators ``Ls`` that generate the liouvillian ``L``.

    A Liouvillian defined by a hermitian Hamilton operator :math:`H` and a
    vector of collapse operators
    :math:`\mathbf{L} = (L_1, L_2, \dots L_n)^T` is invariant under the
    following two operations:

    .. math::
        \left(H, \mathbf{L}\right) & \mapsto \left(H + {1\over 2i}\left(\mathbf{w}^\dagger \mathbf{L} - \mathbf{L}^\dagger \mathbf{w}\right), \mathbf{L} + \mathbf{w} \right) \\
        \left(H, \mathbf{L}\right) & \mapsto \left(H, \mathbf{U}\mathbf{L}\right)\\

    where :math:`\mathbf{w}` is just a vector of complex numbers and
    :math:`\mathbf{U}` is a complex unitary matrix.  It turns out that for
    quantum optical circuit models the set of collapse operators is often
    linearly dependent.  This routine tries to find a representation of the
    Liouvillian in terms of a Hamilton operator ``H`` with as few non-zero
    collapse operators ``Ls`` as possible.  Consider the following example,
    which results from a two-port linear cavity with a coherent input into the
    first port:

    >>> kappa_1, kappa_2 = sympy.symbols('kappa_1, kappa_2', positive = True)
    >>> Delta = sympy.symbols('Delta', real = True)
    >>> alpha = sympy.symbols('alpha')
    >>> H = (Delta * Create(hs=1) * Destroy(hs=1) +
    ...      (sqrt(kappa_1) / (2 * I)) *
    ...      (alpha * Create(hs=1) - alpha.conjugate() * Destroy(hs=1)))
    >>> Ls = [sqrt(kappa_1) * Destroy(hs=1) + alpha,
    ...       sqrt(kappa_2) * Destroy(hs=1)]
    >>> LL = liouvillian(H, Ls)
    >>> Hnf, Lsnf = liouvillian_normal_form(LL)
    >>> print(ascii(Hnf))
    -I*alpha*sqrt(kappa_1) * a^(1)H + I*sqrt(kappa_1)*conjugate(alpha) * a^(1) + Delta * a^(1)H * a^(1)
    >>> len(Lsnf)
    1
    >>> print(ascii(Lsnf[0]))
    sqrt(kappa_1 + kappa_2) * a^(1)

    In terms of the ensemble dynamics this final system is equivalent.
    Note that this function will only work for proper Liouvillians.


    Args:
        L (SuperOperator): The Liouvillian

    Returns:
        tuple: ``(H, Ls)``

    Raises:
        .BadLiouvillianError

def liouvillian_normal_form(L, symbolic = False):
    r"""Return a Hamilton operator ``H`` and a minimal list of collapse
    operators ``Ls`` that generate the liouvillian ``L``.

    A Liouvillian defined by a hermitian Hamilton operator :math:`H` and a
    vector of collapse operators
    :math:`\mathbf{L} = (L_1, L_2, \dots L_n)^T` is invariant under the
    following two operations:

    .. math::
        \left(H, \mathbf{L}\right) & \mapsto \left(H + {1\over 2i}\left(\mathbf{w}^\dagger \mathbf{L} - \mathbf{L}^\dagger \mathbf{w}\right), \mathbf{L} + \mathbf{w} \right) \\
        \left(H, \mathbf{L}\right) & \mapsto \left(H, \mathbf{U}\mathbf{L}\right)\\

    where :math:`\mathbf{w}` is just a vector of complex numbers and
    :math:`\mathbf{U}` is a complex unitary matrix.  It turns out that for
    quantum optical circuit models the set of collapse operators is often
    linearly dependent.  This routine tries to find a representation of the
    Liouvillian in terms of a Hamilton operator ``H`` with as few non-zero
    collapse operators ``Ls`` as possible.  Consider the following example,
    which results from a two-port linear cavity with a coherent input into the
    first port:

    >>> kappa_1, kappa_2 = sympy.symbols('kappa_1, kappa_2', positive = True)
    >>> Delta = sympy.symbols('Delta', real = True)
    >>> alpha = sympy.symbols('alpha')
    >>> H = (Delta * Create(hs=1) * Destroy(hs=1) +
    ...      (sqrt(kappa_1) / (2 * I)) *
    ...      (alpha * Create(hs=1) - alpha.conjugate() * Destroy(hs=1)))
    >>> Ls = [sqrt(kappa_1) * Destroy(hs=1) + alpha,
    ...       sqrt(kappa_2) * Destroy(hs=1)]
    >>> LL = liouvillian(H, Ls)
    >>> Hnf, Lsnf = liouvillian_normal_form(LL)
    >>> print(ascii(Hnf))
    -I*alpha*sqrt(kappa_1) * a^(1)H + I*sqrt(kappa_1)*conjugate(alpha) * a^(1) + Delta * a^(1)H * a^(1)
    >>> len(Lsnf)
    1
    >>> print(ascii(Lsnf[0]))
    sqrt(kappa_1 + kappa_2) * a^(1)

    In terms of the ensemble dynamics this final system is equivalent.
    Note that this function will only work for proper Liouvillians.


    Args:
        L (SuperOperator): The Liouvillian

    Returns:
        tuple: ``(H, Ls)``

    Raises:
        .BadLiouvillianError
    """
    L = L.expand()

    if isinstance(L, SuperOperatorPlus):
        spres = []
        sposts = []
        collapse_form = defaultdict(lambda: defaultdict(int))
        for s in L.operands:
            if isinstance(s, ScalarTimesSuperOperator):
                coeff, term = s.operands
            else:
                coeff, term = One, s
            if isinstance(term, SPre):
                spres.append(coeff * term.operands[0])
            elif isinstance(term, SPost):
                sposts.append((coeff * term.operands[0]))
            else:
                if (not isinstance(term, SuperOperatorTimes) or not
                        len(term.operands) == 2 or not
                        (isinstance(term.operands[0], SPre) and
                        isinstance(term.operands[1], SPost))):
                    raise BadLiouvillianError(
                            "All terms of the Liouvillian need to be of form "
                            "SPre(X), SPost(X) or SPre(X)*SPost(X): This term "
                            "is in violation {!s}".format(term))
                spreL, spostL = term.operands
                Li, Ljd = spreL.operands[0], spostL.operands[0]

                try:
                    complex(coeff)
                except (ValueError, TypeError):
                    symbolic = True
                    coeff = coeff.simplify_scalar()

                collapse_form[Li][Ljd] = coeff

        basis = sorted(collapse_form.keys())

        warn_msg = ("Warning: the Liouvillian is probably malformed: "
                    "The coefficients of SPre({!s})*SPost({!s}) and "
                    "SPre({!s})*SPost({!s}) should be complex conjugates "
                    "of each other")
        for ii, Li in enumerate(basis):
            for Lj in basis[ii:]:
                cij = collapse_form[Li][Lj.adjoint()]
                cji = collapse_form[Lj][Li.adjoint()]
                if cij !=0 or cji !=0:
                    diff = (cij.conjugate() - cji)
                    try:
                        diff = complex(diff)
                        if abs(diff) > 1e-6:
                            print(warn_msg.format(Li, Lj.adjoint(), Lj,
                                                  Li.adjoint()))
                    except ValueError:
                        symbolic = True
                        if diff.simplify():
                            print("Warning: the Liouvillian my be malformed, "
                                  "convert to numerical representation")
        final_Lis = []
        if symbolic:
            if len(basis) == 1:
                l1 = basis[0]
                kappa1 = collapse_form[l1][l1.adjoint()]
                final_Lis = [sqrt(kappa1) * l1]
                sdiff = (l1.adjoint() * l1 * kappa1 / 2)
                spres.append(sdiff)
                sposts.append(sdiff)
#            elif len(basis) == 2:
#                l1, l2 = basis
#                kappa_1 = collapse_form[l1][l1.adjoint()]
#                kappa_2 = collapse_form[l2][l2.adjoint()]
#                kappa_12 = collapse_form[l1][l2.adjoint()]
#                kappa_21 = collapse_form[l2][l1.adjoint()]
##                assert (kappa_12.conjugate() - kappa_21) == 0
            else:
                M = SympyMatrix(len(basis), len(basis),
                                lambda i,j: collapse_form[basis[i]][basis[j]
                                            .adjoint()])

                # First check if M is already diagonal (sympy does not handle
                # this well, for some reason)
                diag = True
                for i in range(len(basis)):
                    for j in range(i):
                        if M[i,j].apply_rules() != 0 or M[j, i].apply_rules != 0:
                            diag = False
                            break
                    if not diag:
                        break
                if diag:
                    for bj in basis:
                        final_Lis.append(
                                bj * sqrt(collapse_form[bj][bj.adjoint()]))
                        sdiff = (bj.adjoint() * bj *
                                 collapse_form[bj][bj.adjoint()]/2)
                        spres.append(sdiff)
                        sposts.append(sdiff)

                # Try sympy algo
                else:
                    try:
                        data = M.eigenvects()

                        for evalue, multiplicity, ebasis in data:
                            if not evalue:
                                continue
                            for b in ebasis:
                                new_L = (sqrt(evalue) * sum(cj[0] * Lj
                                         for (cj, Lj)
                                         in zip(b.tolist(), basis))).expand()
                                final_Lis.append(new_L)
                                sdiff = (new_L.adjoint() * new_L / 2).expand()
                                spres.append(sdiff)
                                sposts.append(sdiff)

                    except NotImplementedError:
                        raise CannotSymbolicallyDiagonalize((
                            "The matrix {} is too hard to diagonalize "
                            "symbolically. Please try converting to fully "
                            "numerical representation.").format(M))
        else:
            M = np_array([[complex(collapse_form[Li][Lj.adjoint()])
                           for Lj in basis] for Li in basis])

            vals, vecs = eigh(M)
            for sv, vec in zip(np_sqrt(vals), vecs.transpose()):
                new_L = sum((sv * ci) * Li for (ci, Li) in zip(vec, basis))
                final_Lis.append(new_L)
                sdiff = (.5 * new_L.adjoint()*new_L).expand()
                spres.append(sdiff)
                sposts.append(sdiff)

        miHspre = sum(spres)
        iHspost = sum(sposts)

        if ((not (miHspre + iHspost) is ZeroOperator) or not
                (miHspre.adjoint() + miHspre) is ZeroOperator):
            print("Warning, potentially malformed Liouvillian {!s}".format(L))

        final_H = (I*miHspre).expand()
        return final_H, final_Lis

    else:
        if L is ZeroSuperOperator:
            return ZeroOperator, []

        raise BadLiouvillianError(str(L))