Returns log(exp(lx) - exp(ly)) without leaving log space.

def log_difference(lx, ly):
    """Returns log(exp(lx) - exp(ly)) without leaving log space."""

    # Negative log of double-precision infinity
    li = -709.78271289338397
    diff = ly - lx
    # Make sure log-difference can succeed
    if np.any(diff >= 0):
        raise ValueError(
            'Cannot compute log(x-y), because y>=x for some elements.')
    # Otherwise evaluate log-difference
    return lx + np.log(1. - np.exp(diff))