function lcs = logsumexp(X, varargin)
% LOGSUMEXP(X, dim) computes log(sum(exp(X), dim)) robustly. Care lightspeed users!
%
%     lse = logsumexp(X[, dim]);
%
% This routine works with general ND-arrays and matches Matlab's default
% behavior for sum: if dim is omitted it sums over the first non-singleton
% dimension.
%
% Note: Tom Minka's lightspeed has a logsumexp function, which:
%     1) sets dim=1 if dim is missing
%     2) returns Inf for sums containing Infs and NaNs;

% Iain Murray, September 2010

% History: IM wrote a bad logsumexp in ~2002, then used Tom Minka's version for
% years until eventually wanting something slightly different.

if (numel(varargin) > 1)
    error('Too many arguments')
end

if isempty(X)
    % Easiest way to get this trivial but annoying case right!
    lcs = log(sum(exp(X),varargin{:}));
    return;
end

if isempty(varargin)
    mx = max(X);
else
    mx = max(X, [], varargin{:});
end
Xshift = bsxfun(@minus, X, mx);

% There is potentially a speed benefit from ignoring terms in the sum that are
% much smaller than the biggest one. But this is problem specific. It's also
% tricky to neatly omit the zero terms from the sum for general arrays and
% summing directions. I decided not to try this optimization in the default
% version.
expXshift = zeros(size(Xshift));
msk = Xshift > -40;
expXshift(msk) = exp(Xshift(msk));

lcs = bsxfun(@plus, log(sum(expXshift,varargin{:})), mx);

idx = isinf(mx);
lcs(idx) = mx(idx);
lcs(any(isnan(X),varargin{:})) = NaN;