def kullback_leibler(h1, h2): # 83 us @array, 109 us @list \w 100 bins
    r"""
    Kullback-Leibler divergence.
    
    Compute how inefficient it would to be code one histogram into another.
    Actually computes :math:`\frac{d_{KL}(h1, h2) + d_{KL}(h2, h1)}{2}` to achieve symmetry.
    
    The Kullback-Leibler divergence between two histograms :math:`H` and :math:`H'` of size
    :math:`m` is defined as:
    
    .. math::
    
        d_{KL}(H, H') = \sum_{m=1}^M H_m\log\frac{H_m}{H'_m}
    
    *Attributes:*

    - quasimetric (but made symetric)
    
    *Attributes for normalized histograms:*

    - :math:`d(H, H')\in[0, \infty)`
    - :math:`d(H, H) = 0`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-normalized histograms:*

    - not applicable
    
    *Attributes for not-equal histograms:*

    - not applicable
        
    Parameters
    ----------
    h1 : sequence
        The first histogram, where h1[i] > 0 for any i such that h2[i] > 0, normalized.
    h2 : sequence
        The second histogram, where h2[i] > 0 for any i such that h1[i] > 0, normalized, same bins as ``h1``.
    
    Returns
    -------
    kullback_leibler : float
        Kullback-Leibler divergence.

    """
    old_err_state = scipy.seterr(divide='raise')
    try:
        h1, h2 = __prepare_histogram(h1, h2)
        result = (__kullback_leibler(h1, h2) + __kullback_leibler(h2, h1)) / 2.
        scipy.seterr(**old_err_state)
        return result
    except FloatingPointError:
        scipy.seterr(**old_err_state)
        raise ValueError('h1 can only contain zero values where h2 also contains zero values and vice-versa')