import numpy as np

############################################################################
# Functions to compute in log-domain.
############################################################################

def logzero():
    return -np.inf

def safe_log(x):
    if x == 0:
        return logzero()
    return np.log(x)

def logsum_pair(logx, logy):
    '''
    Return log(x+y), avoiding arithmetic underflow/overflow.

    logx: log(x)
    logy: log(y)

    Rationale:

    x + y    = e^logx + e^logy
             = e^logx (1 + e^(logy-logx))
    log(x+y) = logx + log(1 + e^(logy-logx)) (1)

    Likewise,
    log(x+y) = logy + log(1 + e^(logx-logy)) (2)

    The computation of the exponential overflows earlier and is less precise
    for big values than for small values. Due to the presence of logy-logx
    (resp. logx-logy), (1) is preferred when logx > logy and (2) is preferred
    otherwise.
    '''
    if logx == logzero():
        return logy
    elif logx > logy:
        return logx + np.log1p(np.exp(logy-logx))
    else:
        return logy + np.log1p(np.exp(logx-logy))

def logsum(logv):
    '''
    Return log(v[0]+v[1]+...), avoiding arithmetic underflow/overflow.
    '''
    res = logzero()
    for val in logv:
        res = logsum_pair(res, val)
    return res


############################################################################
# This implementation is faster, but may give problems with log(0), so I
# commented it out
# def logsum(logv):
#     '''
#     Return log(v[0]+v[1]+...), avoiding arithmetic underflow/overflow.
#     '''
#     c = np.max(logv)
#     return c + np.log(np.sum(np.exp(logv - c)))
############################################################################
