def band_gap(bands, occ, efermi=None):
    """
    Find the band gap in a bandstructure.

    :param numpy.array bands:
        2D bands array (as from BandsData.get_bands())
    :param numpy.array occ:
        2D occupations array matching bands
    :param float efermi:
        optional: Fermi Energy level
    :return:
        {'gap': gap energy difference,
         'direct': bool wether it's a direct band gap
         'vector': start and end points to draw the
         band gap as an arrow in the bandstructure plot.
         the points consist of (k, Dispersion) coordinates
         }
    """
    assert bands.shape == occ.shape
    result = {'gap': None, 'direct': None, 'vector': []}
    nbands = bands.shape[1]
    occupied = [i for i in range(nbands) if occ[:, i].any()]
    unoccupied = [i for i in range(nbands) if not occ[:, i].any()]
    # if either homo or lumo is not included, no info can be given
    if not occupied or not unoccupied:  # if any of them is empty
        return result
    # highest band with any occupation
    homo = bands[:, max(occupied)]
    # lowest completely unoccupied band
    lumo = bands[:, min(unoccupied)]
    gap_lower = homo.max()
    gap_lower_k = homo.argmax()
    # if homo crosses efermi, there is no band gap
    if efermi:
        if gap_lower > efermi:
            result['gap'] = 0
            result['direct'] = False
            return result
    gap_upper = lumo.min()
    gap_upper_k = lumo.argmin()
    # if lumo crosses efermi, something is wrong
    if efermi:
        if gap_upper < efermi:
            raise ValueError(('The given E_fermi was higher than ' 'the lowest point of the lowest unoccupied band'))
    result['gap'] = gap_upper - gap_lower
    result['direct'] = bool(gap_upper_k == gap_lower_k)
    # check wether the two closest points are at the same kpoint (direct)
    # or not (indirect)
    result['vector'] = [(gap_lower_k, gap_lower), (gap_upper_k, gap_upper)]
    return result