# Models log(sum(exp(Ax+b))) <= 0.
# Each row of [A b] describes one of the exp-terms
def logsumexp(M, A, x, b):    
    u = M.variable(A.shape[0])
    M.constraint( Expr.sum(u), Domain.lessThan(1.0))
    M.constraint( Expr.hstack(u,
                              Expr.constTerm(A.shape[0], 1.0),
                              Expr.add(Expr.mul(A, x), b)), Domain.inPExpCone())