docs: def define_simplification(omega_level, xi, Nl):
      """Return a simplifying function, its inverse, and simplified frequencies.
  
      This implements an index iu that labels energies in a non-degenerate
      way.
  
      >>> Ne = 6
      >>> Nl = 2
      >>> omega_level = [0.0, 100.0, 100.0, 200.0, 200.0, 300.0]
      >>> xi = np.zeros((Nl, Ne, Ne))
      >>> coup = [[(1, 0), (2, 0)], [(3, 0), (4, 0), (5, 0)]]
      >>> for l in range(Nl):
      ...     for pair in coup[l]:
      ...         xi[l, pair[0], pair[1]] = 1.0
      ...         xi[l, pair[1], pair[0]] = 1.0
  
      >>> aux = define_simplification(omega_level, xi, Nl)
      >>> u, invu, omega_levelu, Neu, xiu = aux
      >>> print(omega_levelu)
      [0.0, 100.0, 200.0, 300.0]
      >>> print(Neu)
      4
      >>> print(xiu)
      [[[0. 1. 0. 0.]
        [1. 0. 0. 0.]
        [0. 0. 0. 0.]
        [0. 0. 0. 0.]]
      <BLANKLINE>
       [[0. 0. 1. 1.]
        [0. 0. 0. 0.]
        [1. 0. 0. 0.]
        [1. 0. 0. 0.]]]
  
      """
      try:
          Ne = len(omega_level)
      except:
          Ne = omega_level.shape[0]
      #####################################
      # 1 We calculate the symplifying functions.
      om = omega_level[0]
      iu = 0; Neu = 1
      omega_levelu = [om]
      d = {}; di = {0: 0}
      for i in range(Ne):
          if omega_level[i] != om:
              iu += 1
              om = omega_level[i]
              Neu += 1
              omega_levelu += [om]
              di.update({iu: i})
          d.update({i: iu})
  
      def u(i):
          return d[i]
  
      def invu(iu):
          return di[iu]
  
      #####################################
      # 2 We build the simplified xi.
      Neu = len(omega_levelu)
      xiu = np.array([[[xi[l, invu(i), invu(j)]
                      for j in range(Neu)] for i in range(Neu)]
                      for l in range(Nl)])
      #####################################
      return u, invu, omega_levelu, Neu, xiu