function [Am,Vm,bm,Bm,taum,loglikm,mhjumps,logpm] = gibbsvoigt(X,K,M,varargin)
% GIBBSVOIGT Pseudo-voigt source seperation Gibbs sampler
%
% Usage
% [Am,Vm,bm,Bm,taum,loglikm,mhjumps,logpm] = gibbsvoigt(X,N,M,[options])
%
% Input
% X Data matrix (W x N)
% K Number of Voigt components
% M Number of samples to compute
% options
% alpha Prior for A (default 0) formula (7)
% muc Value of prior for c (default zeros(K,1)) formula (8)
% tauc Value of prior for c (default zeros(K,1)) formula (8)
% mug Value of prior for gamma (default 10*ones(K,1)) formula (8)
% taug Value of prior for gamma (default 0.5*ones(K,1)) formula (8)
% mub Value of prior for b (default 0.5) formula (9)
% taub Value of prior for b (default 10) formula (9)
% muB Value of prior for B (default 0.5) formula (9)
% tauB Value of prior for B (default 10) formula (9)
% atau Prior for atau (default 0) formula (12a)
% btau Prior for btau (default 0) formula (12a)
%
% A Initial value for A (K x N)
% V Initial value for V (3 x K)
% b Initial value for A (1 x N)
% B Initial value for B (1 x W)
% tau Initial value for precision (default 1)
%
% chains Number of chains to run (default 1)
% skip Number of initial samples to skip (default 100)
% stride Return every n'th sample (default 1)
% mhstep Steps to use for MH sampler (1 x 3)
%
% gpu Use GPU to calculate inference (logical)
% gpuprec Floating point precision to use on gpu (matlab type string)
%
% nA Do not sample from these columns of A (K x 1 logical)
% nV Do not sample from these voigts V (K x 1 logical)
% nb Do not sample from b (logical)
% nB Do not sample from B (logical)
% ns Do not sample from tau (logical)
%
% Output
% Am Samples of A (K x N x M)
% Vm Samples of A (3 x K x M)
% bm Samples of n (M x N)
% Bm Samples of B (M x W)
% taum Samples of taum (M x 1)
%
% Authors
% Mikkel N. Schmidt,
% Tommy S. Alstrøm
% DTU Compute, Technical University of Denmark.
%
% Method described in the paper
% A pseudo-Voigt component model for high-resolution recovery of
% constituent spectra in Raman spectroscopy, Alstrøm et al, IEEE 2017
% International Conference on Acoustics, Speech, and Signal Processing
% (ICASSP 2017)
% Copyright DTU Compute
[W,N] = size(X);
opts = mgetopt(varargin);
alpha = updim(mgetopt(opts, 'alpha', zeros(K,N)),[K N]);
atau = mgetopt(opts, 'atau', 0);
btau = mgetopt(opts, 'btau', 0);
taub = mgetopt(opts, 'taub', 10);
mub = mgetopt(opts, 'mub', 0.5);
tauB = mgetopt(opts, 'tauB', 10);
muB = mgetopt(opts, 'muB', 0.5);
muc = mgetopt(opts, 'muc', zeros(K,1));
tauc = mgetopt(opts, 'tauc', zeros(K,1));
mug = mgetopt(opts, 'mug', 10*ones(K,1));
taug = mgetopt(opts, 'taug', 0.5*ones(K,1));
A0 = mgetopt(opts, 'A', rand(K,N));
b0 = mgetopt(opts, 'b', rand(1,N));
B0 = mgetopt(opts, 'B', rand(1,W));
V0 = mgetopt(opts, 'V', [W*rand(K,1) repmat([10 0.5],K,1)]);
tau0 = mgetopt(opts, 'tau', 1);
chains = mgetopt(opts, 'chains', 1);
stride = mgetopt(opts, 'stride', 1);
skip = mgetopt(opts, 'skip', 0);
nA = mgetopt(opts, 'nA', false(K,1));
nb = mgetopt(opts, 'nb', false);
nB = mgetopt(opts, 'nB', false);
nV = mgetopt(opts, 'nV', false(K,1));
nt = mgetopt(opts, 'nt', false);
gpuprec = mgetopt(opts, 'gpuprec', 'single');
gpu = mgetopt(opts, 'gpu', false);
mhstep = mgetopt(opts, 'mhstep', [0.1 0.01 0.001]);
Am = zeros(K,N,M*chains);
bm = zeros(M*chains,N);
Bm = zeros(M*chains,W);
Vm = zeros(K,3,M*chains);
logpm = zeros(M*chains,1);
loglikm = zeros(M*chains,1);
taum = zeros(M*chains,1);
WW = 1:W;
mhjumps(K,1) = 0;
muv = [muc mug];
tauv = [tauc taug];
if gpu
Am = gpuArray(cast(Am,gpuprec));
Bm = gpuArray(cast(Bm,gpuprec));
taum = gpuArray(cast(taum,gpuprec));
A0 = gpuArray(cast(A0,gpuprec));
B0 = gpuArray(cast(B0,gpuprec));
X = gpuArray(cast(X,gpuprec));
end
fVoigt = any(~nV);
floglik = fVoigt || ~nt || true;
fXAV = floglik || ~nB || ~nb;
for r = 1:chains
A = A0; % Voigt weigths
b = b0; % Background signal weights
B = B0; % Background signal
V = V0; % Voigt parameters
tau = tau0;
loglik_old = -inf;
loglik = NaN;
VW = voigt(V,WW); % Full voigts
VV = VW*VW';
VX = VW*X;
Bb = B'*b;
bl_ = zeros(1,N);
Al_ = zeros(1,N);
Bl_ = zeros(1,W);
for m = 1:M
for i = 1:skip*(m==1)+stride*(m>1)
VBb = VW*Bb;
for k = 1:K
if ~nA(k)
% sample a's (voigt weights with gibbs)
kk = [1:k-1 k+1:K];
tau_ = tau*VV(k,k);
if K==1
mu_ = (tau*(VX-VBb)-alpha)/tau_;
else
mu_ = (tau*(VX(k,:)-VBb(k,:)-VV(k,kk)*A(kk,:))-alpha(k,:))/tau_;
end
A(k,:) = randr(mu_, 1/tau_, Al_);
end
end
if fXAV
XAV = X-VW'*A;
end
% sample background
if ~nb
tau_ = taub+tau*(B*B');
mu_ = (tau*B*XAV+taub*mub)/tau_;
b = randr(mu_, 1/tau_, bl_);
end
if ~nB
tau_ = tauB+tau*(b*b');
mu_ = (tau*b*XAV'+tauB*muB)/tau_;
B = randr(mu_, 1/tau_, Bl_);
end
if( ~nB || ~nb )
% update background
Bb = B'*b;
end
% sample noise
if floglik
innerloglik = 0.5*sum(sum((XAV-Bb).^2));
loglik = W*N/2*log(tau)-tau*innerloglik;
end
if ~nt
aa = (W*N)/2+atau;
bb = btau + innerloglik;
tau = gamrnd(aa, 1/bb);
end
if fVoigt
% update loglik
loglik_old = -tau*innerloglik;
mhrands = rand(K,1);
end
for k = 1:K
if ~nV(k)
% sample max 10 times before valid sample is obtained
for ii=1:10
Vnew = V(k,:)+randn(1,3).*mhstep;
if Vnew(3)>1 || Vnew(1)>W || any(Vnew<0)
continue;
end
break
end
if Vnew(3)>1 || Vnew(1)>W || any(Vnew<0)
continue;
end
VW_ = VW;
VW_(k,:) = voigt(Vnew,WW); % Full voigts
% prior difference term
prior = sum(-tauv(k,:).*(Vnew(1:2).^2-V(k,1:2).^2 ...
-2*muv(k,:).*(Vnew(1:2)-V(k,1:2))));
loglik_new = -tau*0.5*sum(sum((X-VW_'*A-Bb).^2));
loglik_old = -tau*0.5*sum(sum((X-VW'*A-Bb).^2));
if( exp(loglik_new-loglik_old+prior) > mhrands(k) )
V(k,:) = Vnew;
VW(k,:) = VW_(k,:);
loglik_old = loglik_new;
mhjumps(k) = mhjumps(k) + 1;
else
fail.ratio = exp(loglik_new-loglik_old+prior);
fail.lognew = loglik_new;
fail.logold = loglik_old;
fail.prior = prior;
fail.Vnew = Vnew;
fail.Vold = V(k,:);
end
end
end
end
iRM = m+(r-1)*M;
Am(:,:,iRM) = A;
bm(iRM,:) = b;
Bm(iRM,:) = B;
Vm(:,:,iRM) = V;
taum(iRM) = tau;
logpm(iRM) = loglik_old;
loglikm(iRM) = loglik;
end
end
%--------------------------------------------------------------------------
function X = updim(x, dim)
% UPDIM Replicate and tile an array
%
% Usage
% X = updim(x, dim)
% Copyright 2007 Mikkel N. Schmidt, ms@it.dk, www.mikkelschmidt.dk
dimx = zeros(size(dim));
for k = 1:length(dim), dimx(k) = size(x,k); end
dim(dim==dimx) = 1;
X = repmat(x,dim);
%--------------------------------------------------------------------------
function x = randr(m, s, l)
% RANDR Random numbers from
% p(x)=K*exp(-(x-m)^2/s-l'x), x>=0
%
% Usage
% x = randr(m,s,l)
% Copyright 2007 Mikkel N. Schmidt, ms@it.dk, www.mikkelschmidt.dk
if isa(m,'gpuArray')
m = gather(m);
s = gather(s);
end
A = (l.*s-m)./(sqrt(2*s));
a = A>26;
%else
x = zeros(size(m),'like',m);
y = rand(size(m),'like',m);
%end
x(a) = -log(y(a))./((l(a).*s-m(a))./s);
R = erfc(abs(A(~a)));
x(~a) = erfcinv(y(~a).*R-(A(~a)<0).*(2*y(~a)+R-2)).*sqrt(2*s)+m(~a)-l(~a).*s;
x(isnan(x)) = 0;
x(x<0) = 0;
x(isinf(x)) = 0;
x = real(x);
%--------------------------------------------------------------------------
function out = mgetopt(varargin)
% MGETOPT Parser for optional arguments
%
% Usage
% Get alpha parameter structure from 'varargin'
% opts = mgetopt(varargin);
%
% Get and parse alpha parameter:
% var = mgetopt(opts, varname, default);
% opts: parameter structure
% varname: name of variable
% default: default value if variable is not set
%
% var = mgetopt(opts, varname, default, command, argument);
% command, argument:
% String in set:
% 'instrset', {'str1', 'str2', ... }
%
% Example
% function y = myfun(x, varargin)
% ...
% opts = mgetopt(varargin);
% parm1 = mgetopt(opts, 'parm1', 0)
% ...
% Copyright 2007 Mikkel N. Schmidt, ms@it.dk, www.mikkelschmidt.dk
if nargin==1
if isempty(varargin{1})
out = struct;
elseif isstruct(varargin{1})
out = varargin{1}{:};
elseif isstruct(varargin{1}{1})
out = varargin{1}{1};
else
out = cell2struct(varargin{1}(2:2:end),varargin{1}(1:2:end),2);
end
elseif nargin>=3
opts = varargin{1};
varname = varargin{2};
default = varargin{3};
validation = varargin(4:end);
if isfield(opts, varname)
out = opts.(varname);
else
out = default;
end
for narg = 1:2:length(validation)
cmd = validation{narg};
arg = validation{narg+1};
switch cmd
case 'instrset',
if ~any(strcmp(arg, out))
fprintf(['Wrong argument %sigma = ''%sigma'' - ', ...
'Using default : %sigma = ''%sigma''\n'], ...
varname, out, varname, default);
out = default;
end
case 'dim'
if ~all(size(out)==arg)
fprintf(['Wrong argument dimension: %sigma - ', ...
'Using default.\n'], ...
varname);
out = default;
end
otherwise,
error('Wrong option: %sigma.', cmd);
end
end
end