function y = log1pexp(x)
%LOG1PEXP computes log(1+exp(x)) without under or overflow (USE WITH CARE)
%
% Some aspects of the numerics of this routine aren't well thought out. As with
% all routines, users are encouraged to verify for themselves that log1pexp() is
% suitable for their purposes.

% Iain Murray, November 2007

% Masks for three different methods:
nan_mask = isnan(x);
small = (x < -10); % Slightly arbitrary
naive_mask = ~(nan_mask | small);

% x is fairly big. Just avoid hitting Infs:
% Numerically log(1+exp(x)) == x for large x.
y = zeros(size(x));
y(naive_mask) = log(1 + exp(min(x(naive_mask), 500)));
y(naive_mask) = max(x(naive_mask), y(naive_mask));

% x is fairly small. Use the trick from:
% http://docs.sun.com/source/806-3568/ncg_goldberg.html
% But I'm just justifying the 'small' range empirically. I haven't bothered to
% work out the exact properties of this and can imagine getting burnt in
% the future. Be careful!
simple_exp = exp(x(small));
simple_exp_p1 = simple_exp + 1;
careful = simple_exp.*log(simple_exp_p1)./(simple_exp_p1-1);
mask1 = simple_exp_p1==1;
careful(mask1) = simple_exp(mask1);
y(small) = careful;

y(nan_mask) = NaN;