%=====================================================================
% Programmers: 
% Fa-Yu Wang, E-mail: fayuwang@gmail.com
% Tsung-Han Chan, E-mail: chantsunghan@gmail.com
% Date: 04/27/2009
% -------------------------------------------------------
% Reference: 
% F.-Y. Wang, C.-Y. Chi, T.-H. Chan, and Y. Wang, ¡§Non-negative
% least-correlated component analysis for separation of dependent sources 
% by volume maximization,¡¨to appear in IEEE Trans. Pattern Analysis and 
% Machine Intelligence, 2009 (regular paper).
%==========================================================================
% function [Y iter_cnt time]=nLCA_IVM(X, nor_flag, show_flag)
%--------------------------------------
% Input: 
%--------------------
% X  is the N by L observation matrix, where L is the data length and N is
% the number of sources
% nor_flag: 1- normalize each row of X such that sum(each row of X) equals
% unity, and 0 - do nothing. 
% show_flag: 1- display current information in nLCA-IVM, and 0 - otherwise 
%==========================================================================
% Output: 
%--------------------
% Y is the N by L extracted source matrix
% iter_cnt is the total number of iterations passed by nLCA-IVM
% time: computation time (seconds)
%==========================================================================

function [Y iter_cnt time]=nLCA_IVM(X, nor_flag, show_flag)

t0 = clock;
%--------Row sum-based normalization to X--------------
if nor_flag==1
    X=X./(sum(X')'*ones(1,size(X,2)));
end

X_ori = X;
N = size(X_ori, 1);
% ----- reduce redundant constraint set by Quick hull algorithm--------------
if N < 9
    X=red_const_qhull(X')'; 
% The procedure may not be available for N > = 9 due to our of memory problem 
% which is a the intrinsic problem in quick hull algorithm. 
else 
    X = X;
end
%-------------
TOL_convergence = 1e-5;  % convergence tolarance 
Iteration_max=100; % the maximum number of iterations
[L]=size(X, 2);
iter_cnt=0;
rec=100; 
%-------------
W=eye(N); % initial unmixing matrix
while ((iter_cnt<Iteration_max)&(rec>TOL_convergence))
    iter_cnt=iter_cnt+1;
    obj0=abs(det(W));
    for row=1:N
        other=find(row~=1:N);
        temp=W(other,:);
        for k=1:N
             cc(k)=(-1)^(row+k)*det(temp(:,find(k~=1:N)));
        end  
        %--relax the absolutely objective function and add an extra constraint-------
        %%%%1%%%%
        B1=sparse(-[X';cc]);
        c=sparse([cc';zeros(L+1,1)]);
        ttt=sparse(1:L+1,1:L+1,ones(L+1,1),L+1,L+1);
        A1=sparse([B1 ttt;ones(1,N) zeros(1,L+1)]);
        b=sparse([zeros(L+1,1);1]);
    
        K.f=N;
        K.l=L+1;
        pars.fid=0;
        x1=sedumi(A1,b,c,K,pars);
        w1=x1(1:N);
        %%%%2%%%%
        B2=sparse([-X';cc]);
        A2=sparse([B2 ttt;ones(1,N) zeros(1,L+1)]);
        x2=sedumi(A2,b,c,K,pars);
        w2=x2(1:N);
        %----compare---------------------
        w=[w1 w2];
        [val ind]=max([cc*w1 -cc*w2]);
        w_true=w(:,ind);
        %----update W------
        W(row,:) = w_true';
    end
    rec = abs(obj0-abs(det(W)))/obj0;
    if show_flag
        disp(' ');
        disp(strcat('Number of iterations: ', num2str(iter_cnt)))
        disp(strcat('Relative change in objective: ', num2str(rec)))
        disp(strcat('Current volume: ', num2str(abs(det(W)))))
    end
end
Y = W*X_ori;


%----rescale Y such that max{Y}=1
Y=Y./(max(Y')'*ones(1,size(Y,2)));

time = etime(clock,t0);


% ============================================== %
% Quick hull algorithm
% ============================================== %
function Xf=red_const_qhull(X)

[L,M]=size(X);
index=find(sum(abs(X'))>=1e-6);
Xn=X(index,:); 
LL=length(index);
Y=sum(Xn');
Xh=Xn./(Y'*ones(1,M));
K = convhulln(Xh,{'QJ'});
index_s=sort(vec(K)); % sorting
[leng,t]=size(index_s);
p=zeros(leng,1);
p(1,1)=1;
l=1;
while (l<leng)
    if index_s(l+1)~=index_s(l)
       p(l+1,1)=1;
    end
    l=l+1;
end
index_f=index_s(find(p>0)); % list the indices of data corresponding to possible constraints
Xf=Xh(index_f,:);