def log_sum_exp_ab(a, b):
    """log_sum_exp for two scalars.

    Parameters
    ----------
    a, b: float 

    Returns
    -------
    c: float
        c = log(e^a + e^b)
    """
    if a > b:
        return a + np.log(1. + np.exp(b - a))
    else:
        return b + np.log(1. + np.exp(a - b))