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``.

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