% P = potLaplace(s) - Laplace (double exponential) potential
%
% pot(s) = exp(-|s|)
%
%   See also POTFUNCTIONS.M.
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2010 November 03

function P = potLaplace(s,type,z)

if nargin==1
  q = numel(s); P = zeros(q,4);                   % allocate memory, b=0, d2lp=0
  P(:,1) = -abs(s(:));                                           % log potential
  P(:,2) = 1-2*(s(:)>0);                       % 1st derivative of log potential
else
  if strcmp(type,'VB')
    P = potLaplace(s);
  elseif strcmp(type,'EP')
    q = numel(s); P = zeros(q,3);                 % allocate memory, b=0, d2lp=0
    fac = 1e3;          % factor between the widths of the two distributions ...
       % ... from when one considered a delta peak, we use 3 orders of magnitude
    idpot = fac<sqrt(z);                             % Potential is a delta peak
    idgau = sqrt(z)<1/fac;                            % Gaussian is a delta peak
    id = ~idgau & ~idpot;                          % interesting case in between
    if any(idpot)
      Q = potGauss(s(idpot)./sqrt(z(idpot))); P(idpot,:) = Q(:,1:3);
      P(idpot,1) = P(idpot,1) - log(2*pi*z(idpot))/2 + log(2);
      P(idpot,2) = P(idpot,2)./sqrt(z(idpot));
      P(idpot,3) = P(idpot,3)./z(idpot);
    end
    if any(idgau)
      Q = potLaplace(s(idgau)); P(idgau,:) = Q(:,1:3);
    end
    if any(id)
      % substitution to obtain unit variance, zero mean Laplacian
      tmu = s(id)./sqrt(2); tvar = z(id)./2;

      % an implementation based on logphi(t) = log(normcdf(t))
      zp = (tmu+sqrt(2)*tvar)./sqrt(tvar);
      zm = (tmu-sqrt(2)*tvar)./sqrt(tvar);
      ap =  logphi(-zp)+sqrt(2)*tmu;
      am =  logphi( zm)-sqrt(2)*tmu;
      P(id,1) = logsum2exp([ap,am]) + tvar;

      lqp = -zp.^2/2 - log(2*pi)/2 - logphi(-zp);       % log( N(z)/Phi(z) )
      lqm = -zm.^2/2 - log(2*pi)/2 - logphi( zm);
      dap = -exp(lqp-log(z(id))/2) + 1;
      dam =  exp(lqm-log(z(id))/2) - 1;
                        % ( exp(ap).*dap + exp(am).*dam )./( exp(ap) + exp(am) )
      P(id,2) = expABz_expAx([ap,am],[1;1],[dap,dam],[1;1]);
 
      a = 2./sqrt(z(id));
      bp = 1 - (a - zp./z(id)).*exp(lqp);
      bm = 1 - (a + zm./z(id)).*exp(lqm);
      % d2lZ(id) = ( exp(ap).*bp + exp(am).*bm )./( exp(ap) + exp(am) ) ...
      %            - dlZ(id).^2;
      P(id,3) = expABz_expAx([ap,am],[1;1],[bp,bm],[1;1]) - P(id,2).*P(id,2);
    end
  else
    error('Unknown type')
  end
end

% computes y = log( sum(exp(x),2) ) in a numerically safe way by subtracting 
%  the row maximum to avoid cancelation after taking the exp
%  the sum is done along the rows
function [y,x] = logsum2exp(logx)
  N = size(logx,2);
  max_logx = max(logx,[],2);
  % we have all values in the log domain, and want to calculate a sum
  x = exp(logx-max_logx*ones(1,N));
  y = log(sum(x,2)) + max_logx;

function y = expABz_expAx(A,x,B,z)
  N = size(A,2);  maxA = max(A,[],2);      % number of columns, max over columns
  A = A-maxA*ones(1,N);                                 % subtract maximum value
  y = ( (exp(A).*B)*z ) ./ ( exp(A)*x ); 

% Safe computation of logphi(z) = log(normcdf(z)).
% The function is based on index 5725 in Hart et al. and gsl_sf_log_erfc_e.
function lp = logphi(z)
  z = real(z);                                 % support for real arguments only
  lp = zeros(size(z));                                         % allocate memory
  id1 = z.*z<0.0492;                                 % first case: close to zero
  lp0 = -z(id1)/sqrt(2*pi);
  c = [ 0.00048204; -0.00142906; 0.0013200243174; 0.0009461589032;
       -0.0045563339802; 0.00556964649138; 0.00125993961762116;
       -0.01621575378835404; 0.02629651521057465; -0.001829764677455021;
       2*(1-pi/3); (4-pi)/3; 1; 1];
  f = 0; for i=1:14, f = lp0.*(c(i)+f); end, lp(id1) = -2*f-log(2);
  id2 = z<-11.3137;                                    % second case: very small
  r = [ 1.2753666447299659525; 5.019049726784267463450;
        6.1602098531096305441; 7.409740605964741794425;
        2.9788656263939928886 ];
  q = [ 2.260528520767326969592;  9.3960340162350541504;
       12.048951927855129036034; 17.081440747466004316; 
        9.608965327192787870698;  3.3690752069827527677 ];
  num = 0.5641895835477550741; for i=1:5, num = -z(id2).*num/sqrt(2) + r(i); end
  den = 1.0;                   for i=1:6, den = -z(id2).*den/sqrt(2) + q(i); end
  e = num./den; lp(id2) = log(e/2) - z(id2).^2/2;
  id3 = ~id2 & ~id1; lp(id3) = log(erfc(-z(id3)/sqrt(2))/2);  % third case: rest