function [w_vb, bt]= mpLaplace_truncatedGaussian_vB(Phi, y, Phi_gram, MaxIter, const)
% Function 
% Variational Bayes for linear regression based on a multiparameter 
% truncated Laplace prior. For details see the technical report with title
% “Semisupervised hyperspectral image unmixing using a variational Bayes 
% algorithm”, National Observatory of Athens, June, 2014
% 
% Input
% Phi:      Endmember matrix
% y:        Pixel's spectral measurements
% Phi_gram: Gramian matrix of Phi, Phi.' * Phi 
% MaxIter:  Maximum Iterations ~ 1e3
% const:    small constant ~ 1e-8
% 
% Output
% w_vb:     Estimated abundance vector
% bt:       Estimated noise precision
% 
% Copyright K. Themelis, A. Rontogiannis, K. Koutroumbas

[M,N] = size(Phi);

% complementary indexes
indx = zeros(N-1,N); for n = 1 : N, tmp = 1:N; tmp(n) = []; indx(:,n) = tmp; end

alph = ones(N,1); % initialization
bt = 1;
bn = ones(N,1);
mu = zeros(N,1);

t = 0;
while 1
    t = t + 1;
    mu_old = mu;
    
    % untruncated variance
    sgmi = bt^-1 .* (alph + diag(Phi_gram) ).^-1;
    sgmisd = sqrt(sgmi);
    
    for n = 1 : N
        % untruncated mean
        mu(n) =  Phi(:,n).' * (y - Phi(:,indx(:,n)) * mu(indx(:,n))) / (alph(n) + Phi_gram(n,n));
        
        % truncated mean
        tmp = 1 - .5 * erfc(mu(n) /sqrt(2) /sgmisd(n) );
        if tmp 
            mu(n) = mu(n) + sgmisd(n) / sqrt(2 * pi) * exp(-.5 * mu(n)^2 / sgmi(n) ) / tmp;
        end
    end
    
    % truncated variance
    tmp1 = mu ./ sgmisd; tmp2 = tmp1.^2; tmp3 = 1 - .5 .* erfc(tmp1 ./ sqrt(2)); tmp4 = exp(-tmp2 ./ 2);
    sgmitr = sgmi .* (1 - tmp1 ./ sqrt(2*pi) .* tmp4 ./ tmp3 - (tmp4 ./ sqrt(2*pi) ./ tmp3 ).^2 );
    
    bt = (2e-6 + N + M ) / (2e-6 + sum(alph .* (mu.^2 + sgmitr)) + norm(y-Phi*mu) + sgmitr.' * diag(Phi_gram) );
    
    alph = sqrt( bn ./ (bt * (mu.^2 + sgmitr)));
    
    bn = (2e-6 + 2) ./ (2e-6 + 1./alph + 1./bn );
    
    if (abs(norm(mu - mu_old)) < const || t > MaxIter), w_vb = mu; break; end
end