function S = logsumexp(X,D)
% LOGSUMEXP Log of sum of exponentiated terms
%
% S = LOGSUMEXP(X) is a more numerically stable version of log(sum(exp(X)))
% S = LOGSUMEXP(X,D) takes the sum along dimension D.

% Jerod Weinman jerod@acm.org
%
% $Id: logsumexp.m 45 2012-05-02 17:39:57Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% /home/data/svn/vision/tools/package/GPL_STATEMENT
 
  
  if nargin<2
    D = min(find(size(X)~=1)); % First non-singleton dimension
  end;
  
  if isempty(D)
    D = 1;
  end;

  % Calculate max along dimension
  mx = max(X,[],D);

  % Log of sum of exponentiated data with maximal terms along D
  % subtracted for stability
  S = mx + log( sum( exp( bsxfun(@minus, X, mx)), D));

  S(isinf(mx) & mx<0) = -Inf;
