function [Xk_SL, k_SL, S] = sfa(X)
% sfa.m  significant factor analysis - a program designed
%        to help determine the number of significant factors in a matrix.
%
% Factor Analysis in Chemistry, Third Edition, p387-389
% Edmund R. Malinowki
%
% IND: Indicator Function
% SL: Significance Level
%
% [Xk_SL, k_SL, S] = sfa(X)
% Xk_SL: denoised data matrix with SL values
% k_SL: number of singular values keeped using SL
% S: singular values matrix
% X: data matrix


format short e
max_err = 5;        % probability of being noise for SL

% transpose X if needed
[m,n] = size(X);
transp = 0;
if m < n
    X = X';
    transp = 1;
    [m,n] = size(X);
end

% disp('SVD in progress, please be patient')
% tic
[U,S,V] = svd(X,0);
% toc

%% Indicator function IND
% preallocating
ev = zeros(1,n);
df = zeros(1,n);
rev = zeros(1,n);
sev = zeros(1,n-1);
sdf = zeros(1,n-1);
re = zeros(1,n-1);
ind = zeros(1,n-1);

for j = 1:n
    ev(j) = S(j,j) * S(j,j);
    df(j) = (m-j+1)*(n-j+1);
    rev(j) = ev(j) / df(j);
end

for k = 1:n-1
    sev(k) = sum(ev(k+1:n));
    sdf(k) = sum(df(k+1:n));
end

for i = 1:n-1
    re(i) = sqrt(sev(i) / (m * (n-i)));     % see eq. 4.44
    ind(i) = re(i) / (n-i)^2;               % see eq. 4.63
end

[~,k_IND ] = min(ind);
% disp(['IND function indicates ', int2str(k_IND ), ' significant factors.'])
% disp(['The real error (RE) is +/-', num2str(re(k_IND )), '.'])

%% Significance Level SL
re(n) = NaN; ind(n) = NaN;
t = zeros(n,5);
for j = 1:n
    t(j,1) = j;
    t(j,2) = ev(j);
    t(j,3) = re(j);
    t(j,4) = ind(j);
    t(j,5) = rev(j);
end

for j = 1:n-1
    f = (sdf(j) * ev(j)) / ((m-j+1) * (n-j+1) * sev(j));
    % convert f (see eq. 4.83) into percent significance level
    if j < n
        tt = sqrt(f);
        df = n - j;
        a = tt / sqrt(df);
        b = df / (df + tt * tt);
        im = df - 2;
        jm = df - 2 * fix(df / 2);
        ss = 1;
        cc = 1;
        ks = 2 + jm;
        fk = ks;
        if (im - 2) >= 0
            for k = ks:2:im
                cc = cc * b * (fk-1) / fk;
                ss = ss + cc;
                fk = fk + 2;
            end
        end
        if (df - 1) > 0
            c1 = .5 + (a * b * ss + atan(a)) * .31831;
        else
            c1 = .5 + atan(a) * .31831;
        end
        if jm <= 0
            c1 = .5 + .5 *a * sqrt(b) * ss;
        end
    end
    s1 = 100 * (1 - c1);
    s1 = 2 * s1;
    if s1 < 1e-2, s1 = 0; end
    t(j,6) = s1;
end
t(n,6) = NaN;

k_SL = find(t(:,6) < max_err, 1, 'last');      % see eq. 4.83
if isempty(k_SL), k_SL = 0; end

% disp(['SL function indicates ', int2str(k_SL), ' significant factors.'])
% if k_SL ~= 0
    % disp(['The real error (RE) is +/-', num2str(re(k_SL)), '.'])
% end

% matrix reconstruction
Xk_IND= U(:,1:k_IND)*S(1:k_IND,1:k_IND)*V(:,1:k_IND)';% using IND
Xk_SL = U(:,1:k_SL)*S(1:k_SL,1:k_SL)*V(:,1:k_SL)';    % using SL

% transpose
if transp == 1
    Xk_IND = Xk_IND';
    Xk_SL = Xk_SL';
end
