wolffd@0: function s = logsumexp(a, dim)
wolffd@0: % Returns log(sum(exp(a),dim)) while avoiding numerical underflow.
wolffd@0: % Default is dim = 1 (rows) or dim=2 for a row vector
wolffd@0: % logsumexp(a, 2) will sum across columns instead of rows
wolffd@0: 
wolffd@0: % Written by Tom Minka, modified by Kevin Murphy
wolffd@0: 
wolffd@0: if nargin < 2
wolffd@0:   dim = 1;
wolffd@0:   if ndims(a) <= 2 & size(a,1)==1
wolffd@0:     dim = 2;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % subtract the largest in each column
wolffd@0: [y, i] = max(a,[],dim);
wolffd@0: dims = ones(1,ndims(a));
wolffd@0: dims(dim) = size(a,dim);
wolffd@0: a = a - repmat(y, dims);
wolffd@0: s = y + log(sum(exp(a),dim));
wolffd@0: %i = find(~finite(y));
wolffd@0: %if ~isempty(i)
wolffd@0: %  s(i) = y(i);
wolffd@0: %end