function [varepsilon unstable] = getEstBound(A, W)
[N,M,K]=size(W);
varepsilon=[];
unstable=0;
for k=1:K
    [varTmp unsTmp IdxRows IdxCols] = getEstBoundPerDataSet(A(:,:,k), W(:,:,k));
    varepsilon=[varepsilon varTmp];
    unstable=unstable+unsTmp;
    if k>=2
        unstable=unstable+sum(abs(IdxRows-IdxRowsPre))+sum(abs(IdxCols-IdxColsPre));
    end
    IdxRowsPre=IdxRows;
    IdxColsPre=IdxCols;
end
unstable=unstable>0;

function [varepsilon unstable IdxRows IdxCols] = getEstBoundPerDataSet(A, W)
% return the ISR for NonGaussian i.i.d. case

N=size(A,1);
estI=W*A;
varepsilon=[];
for i=1:N
    [MaxRow IdxMaxRow]=max(abs(estI(i:N,i:N)));
    [MaxCol IdxMaxCol]=max(MaxRow);
    IdxRow=IdxMaxRow(IdxMaxCol)+i-1;
    IdxCol=IdxMaxCol+i-1;
    IdxRows(i)=IdxRow;
    IdxCols(i)=IdxCol;
    estI([i IdxRow],:)=estI([IdxRow i],:);
    estI(:, [i IdxCol],:)=estI(:, [IdxCol i]);
    estI(i,:)=estI(i,:)/estI(i, i);
    varepsilon=[varepsilon estI(i,(estI(i,:)~=estI(i, i)))];
end;
unstable=(sum(sum(abs(estI)>=1))>N);

j=1;

