function y = softmax(a, beta)
%
% y = softmax(a, beta)
% y(i,j) = exp(a(i,j) * beta) / sum_j exp(a(i,j) * beta)  columns contain probs.
% This is also called a multinomial logit function.
% If beta is omitted, it defaults to 1.
% As beta->inf, softmax -> max. As beta -> 0, softmax -> uniform.

if nargin==1, beta = 1; end

x = exp(beta*a);

if any(isinf(x)) % beta is too high
  % create a delta function
  y = zeros(1,length(a));
  y(argmax(a)) = 1;
else
  y = x ./ repmat(sum(x,1), size(x,1), 1);
  %y = normalise(x);
end
