function [B, C, iter, Z, Fit_err] = CSR_ALS(G, A_ini,B_ini,C_ini,L,Maxiter )
% Corrupt Slice Robust Alternating Least Squars (CSR-ALS)
% reference: the paper

% The formulation is
% min_{A B C Z} sum_i^K ||G^{(3)}-KhatriRaoProduct(C,B)A^T-Z||_F^2
% subject to    ||Z||_{col-0}<=L
% the model is that G^{(3)} in KJxI has some corrupted columns
% reference: 
% X. Fu, N. D. Sidiropoulos, J. H. Tranter, and W.-K. Ma, "A factor analysis framework for power spectra separation and multiple emitter localization,"  IEEE Transactions on Signal Processing, accepted, June, 2015


% input:
% initializations of A B C
% A_ini: I x R 
% B_ini: J x R
% C_ini: K x R
% G is the data, which contains some outlying slabs/slices
% the clean slabs should be like G(q,:,:) = B*Dq(A)*C';  % where Dq(A) means a diagonal matrix holding the qth row
% of A
% R: the rank of this tensor
% L: the number of outliers we belive (can be an upper bound)
% Maxiter: maximal iteration


% Programmed by Xiao Fu: xfu@umn.edu
% First version: Jan 27 2014
% Begins at the living room of Dennis' house, comfortabaly


A = A_ini;  [I, R]=size(A);
B = B_ini;  [J, R]=size(B);
C = C_ini;  [K, R]=size(C);

tol = 1e-6;

[~,~,Ztilde]=constructXYZ(G,R,I,J,K);  

% Xtilde IK x J 
% Ytilde IJ x K
% Ztilde KJ x I  % WE ASSUME THAT THE OUTLIERS ARE THE COLUMNS of Ztilde

X_o = Ztilde; 
Z = zeros(K*J,I);
for iter=1:Maxiter
    
    %------------      update Z ---------------------
    Kr = kr_product(B,C)*A.';
    z_err= sum(abs(X_o - Kr).^2);  
    
    fit(iter)=sum(sum(abs(X_o - Z - Kr).^2)); 
    Z = zeros(K*J,I);
    [~,index_err]=sort(z_err,'descend');
    index_delete = index_err(:,1:L);
     
    index_update = [1:1:I];
    Z(:,index_update)=zeros(K*J,I);
    Z(:,index_delete) = X_o(:,index_delete)-Kr(:,index_delete);
    

    % --------- construct a new G using X_o - Z;
    
    index_delete = sort(index_delete ,'ascend');
    Ztilde_new = X_o - Z;
    for ii=1:length(index_delete)
        G_replace{ii} = reshape(Ztilde_new(:,index_delete(ii)),K,J);
    end
    for k=1:K
        G_del = G(:,:,k);
        for ii=1:length(index_delete)
            G_del(index_delete(ii),:)= G_replace{ii}(k,:); % The impact of the original tensor;
        end
        G_update(:,:,k)=G_del;
    end
       
     
    % ---------------- run anthother PARAFAC -------
   [A,C,B]=TALS(G_update,R,A,C,B,1);

   
   
  if iter>1&&abs(fit(iter)-fit(iter-1))/fit(iter-1)<tol;
      break
  end
   
   
   
end 

Fit_err = fit(end);

end   
  
function [Xtilde,Ytilde,Ztilde] = constructXYZ(X,M,K,N,P)

Xtilde = [];
ap1 = zeros(K,N);
for p=1:1:P,
    for k=1:1:K,
        for n=1:1:N,
            ap1(k,n) = X(k,n,p);
        end
    end
    Xtilde = [Xtilde; ap1];
end

Ytilde = [];
ap2 = zeros(N,P);
for k=1:1:K,
    for n=1:1:N,
        for p=1:1:P,
            ap2(n,p) = X(k,n,p);
        end
    end
    Ytilde = [Ytilde; ap2];
end

Ztilde = [];
ap3 = zeros(P,K);
for n=1:1:N,
    for p=1:1:P,
        for k=1:1:K,
            ap3(p,k) = X(k,n,p);
        end
    end
    Ztilde = [Ztilde; ap3];
end

end

function C = kr_product(A,B)

%A and B must have same number of columns
[N,K] = size(A);
[M,K] = size(B);
C = zeros(M*N,K);
for i1 = 1:K
    C(:,i1) = kron_mod(A(:,i1),B(:,i1));
end



end

function X = kron_mod(A,B)
%KRON Kronecker product.
%   kron(A,B) returns the Kronecker product of two matrices A and B, of
%   dimensions I-by-J and K-by-L respectively. The result is an I*K-by-J*L
%   block matrix in which the (i,j)-th block is defined as A(i,j)*B.

%   Version: 06/02/2011
%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)

[I J] = size(A);
[K L] = size(B);

if ~issparse(A) && ~issparse(B)
    
    % Both matrices are dense.
    A = reshape(A,[1 I 1 J]);
    B = reshape(B,[K 1 L 1]);
    X = reshape(bsxfun(@times,A,B),[I*K J*L]);
    
else
    
    % One of the matrices is sparse.
    [ia,ja,sa] = find(A);
    [ib,jb,sb] = find(B);
    ix = bsxfun(@plus,K*(ia(:)-1).',ib(:));
    jx = bsxfun(@plus,L*(ja(:)-1).',jb(:));
    
    % The @and operator is slightly faster for logicals.
    if islogical(sa) && islogical(sb)
        X = sparse(ix,jx,bsxfun(@and,sb(:),sa(:).'),I*K,J*L);
    else
        X = sparse(ix,jx,double(sb(:))*double(sa(:).'),I*K,J*L);
    end
    
end
end