function [ A_est, S_est, time_volmin ] = VolMin_AO( R_whiten,K,B_pre,X) 

%Input:
%R_whiten: Local Covariances (whitened)
%K: number of sources
%L: frame length
%B: prewhiten matrix
%X: data
%output:
% A_est: channel estimation
% S_est: source estimation
% H_est = KR(conj(A_est),A_est)
% Y_normalized: for initializing VolMin-AO


% copyrights: Xiao Fu
% xfu@umn.edu
% See: Xiao Fu, Wing-Kin Ma, Kejun Huang, and Nicholas Sidiropoulos, 
% "Blind Separation of Quasi-Stationary Sources: Exploiting Convex Geometry
% in Covariance Domain" submitted to IEEE Trans. Signal Process., 2014


%----------------------- Initialization by ProSPA ----------------------

[ A_SPA, S_SPA, timeSPA, B, Y] = ProSPA(R_whiten,K,B_pre,X);

tic;
[n,M] = size(Y); 

Max_iter = 500;

for it = 1:Max_iter
    %---------------- Compute the update of Theta --------------------
    
    [Theta ]= SimplexProj( (B'*Y).' );
    Theta = Theta.';
    YD = Y*Theta';
    [U S V]=svd(YD);
    Bt = V*U(:,1:K)';
    B = Bt';
    
    %--------------- For observation -------------------------------
    obj_ob(it) = sum(sum((Y-B*Theta).^2));
    
    %---------------------------Stopping Criterion -------------------------
    if it>=2&& ( abs(obj_ob(it)-obj_ob(it-1))) <1e-2%||(sum(sum(abs(B-B_pre).^2)))<1e-4;%
        break;
    end
    
    
end

H_est = B;

for k=1:K
    A_temp = reshape(H_est(:,k),sqrt(n),sqrt(n));
    [U1 S1] = eig((A_temp+A_temp')/2);
    A_est(:,k) = U1(:,end);
end
A_est = B_pre*A_est;
S_est = pinv(A_est)*X;
time_volmin =  toc;
time_volmin = time_volmin +  timeSPA;
end





function X = SimplexProj(Y) 
% Projection to a simplex; the codes are from the reference [41]
% W. Wang and M. A. Carreira-Perpinán, “Projection onto the probability
% arXiv preprint arXiv:1309.1541v1, 2013.
[N,D] = size(Y);
X = sort(Y,2,'descend');
Xtmp = (cumsum(X,2)-1)*diag(sparse(1./(1:D)));
X = max(bsxfun(@minus,Y,Xtmp(sub2ind([N,D],(1:N)',sum(X>Xtmp,2)))),0);
end

