docs: 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))