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))
def liouvillian_normal_form ( L , symbolic = False ) : 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 ) else : M = SympyMatrix ( len ( basis ) , len ( basis ) , lambda i , j : collapse_form [ basis [ i ] ] [ basis [ j ] . adjoint ( ) ] ) 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 ) 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 ) )
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