docs: def gap(lamda, sigma, mode='individual'):
      """Compute spectral gap.
  
      Useful for eigenvalue/eigenvector bounds. Computes the gap
      :math:`\delta\geq 0` between two sets of real numbers ``lamda`` and
      ``sigma``. The gap can be computed in several ways and may not exist, see
      the ``mode`` parameter.
  
      :param lamda: a non-empty set
        :math:`\Lambda=\{\lambda_1,\ldots,\lambda_n\}` given as a single real
        number or a list or ``numpy.array`` with real numbers.
      :param sigma: a non-empty set :math:`\Sigma=\{\sigma_1,\ldots,\sigma_m\}`.
        See ``lamda``.
      :param mode: (optional). Defines how the gap should be computed. May be one
        of
  
        * ``'individual'`` (default):
          :math:`\delta=\min_{\substack{i\in\{1,\ldots,n\}\\\\j\in\{1,\ldots,m\}}} |\lambda_i - \sigma_j|`.
          With this mode, the gap is always be defined.
        * ``'interval'``: determine the maximal :math:`\delta` such that
          :math:`\Sigma\subset\mathbb{R}\setminus[\min_{\lambda\in\Lambda}\lambda-\delta,\max_{\lambda\in\Lambda}\lambda+\delta]`.
          If the gap does not exists, ``None`` is returned.
  
      :return: :math:`\delta` or ``None``.
      """
      # sanitize input
      if numpy.isscalar(lamda):
          lamda = [lamda]
      lamda = numpy.array(lamda)
      if numpy.isscalar(sigma):
          sigma = [sigma]
      sigma = numpy.array(sigma)
  
      if not numpy.isreal(lamda).all() or not numpy.isreal(sigma).all():
          raise ArgumentError('complex spectra not yet implemented')
  
      if mode == 'individual':
          return numpy.min(numpy.abs(numpy.reshape(lamda, (len(lamda), 1))
                           - numpy.reshape(sigma, (1, len(sigma)))))
      elif mode == 'interval':
          lamda_min, lamda_max = numpy.min(lamda), numpy.max(lamda)
          # determine all values in sigma<lamda_min or >lamda_max
          sigma_lo = sigma <= lamda_min
          sigma_hi = sigma >= lamda_max
          # is a sigma value in lamda interval?
          if not numpy.all(sigma_lo + sigma_hi):
              return None
          delta = numpy.Infinity
          if numpy.any(sigma_lo):
              delta = lamda_min - numpy.max(sigma[sigma_lo])
          if numpy.any(sigma_hi):
              delta = numpy.min([delta, numpy.min(sigma[sigma_hi]) - lamda_max])
  
          return delta