function A_SymEfrb = SymEfrb(X,Wg,iD,nonlin,nboot)
%%%%%%%%%%%%%%%%%%%%%%%%      Algorithm 2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the Efrb method for Symmetric FastICA, proposed in:
% Enhanced bootstrap method for statistical inference in the ICA model 
% S. Basiri, E. Ollila, V. Koivunen, Signal Processing 138 (2017), 53-62.
% Code by Shahab Basiri, Aalto University 2017
% shahab.basiri@aalto.fi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input:
% X     : (d*n) matrix of whitened data.
% Wg    : (d*d) Symmetric FastICA estimate of demixing matrix 
%         for whitened data
% iD    : (p*d) Moore-Penrose generalized inverse of the whitening matrix
% nonlin: FastICA nonlinearity: e.g. {'pow3'}, {'tanh'}, {'gauss'}, ...
% nboot : The number of required bootstrap samples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output:
% A_DefEfrb : (p*d*nboot) Symmetric EFRB estimates of mixing matrix 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[d, n] = size(X);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A_SymEfrb = zeros(d,d,nboot);
%% Centering:
%%%%%%%%%%%%%%
X = bsxfun(@minus, X, mean(X,2));
%% Computing the matrix of partial derivatives:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = Wg*X;
[gS, EdgS, ~] = gdg(S,nonlin);
m = X*gS'./n;
B = EdgS;
Lam = sqrt(m'*m);

%% Generate bootstrap indices:
%%%%%%%%%%%%%%%%%%
    boot = bootrsp(1:n, nboot)';
%% Bootstrap iteration starts here:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for b = 1:nboot
    %% (1)Construct the bootstrap sample
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    X_star = X(:,boot(b,:));
    %% (2)Centering:
    %%%%%%%%%%%%%%
    X_star = bsxfun(@minus, X_star, mean(X_star,2));
    
    % Covariance Matrix:
    %%%%%%%%%%%%%%%%%%%%%
    C_star = X_star*X_star'./n;
    %% (3)Computing the one-step estimator:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    S_star = Wg*X_star;
    [gS_star, ~, ~] = gdg(S_star,nonlin);
    m_star = X_star*gS_star'./n;
    W_star = zeros(d,d);
    for k = 1:d,
        %% (4)Computing the one-step estimator:
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        m_star_k = m_star(:,k);
        lam_star_k = sqrt(m_star_k'/C_star*m_star_k);
        w1_star_k =C_star\m_star_k./lam_star_k;
        %% (5)Computing the EFRB replica:
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        w_k = Wg(k,:)';
        lam_k = Lam(k,k);
        B_k = B(k);        
        wR_star_k = w1_star_k - w_k*(w_k'*w1_star_k*B_k/lam_k);
        W_star(k,:) = wR_star_k';
    end
    %% (6)Symmetric orthogonalization:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    W_star = ((W_star*C_star*W_star')^(-1/2))*W_star;
    %% Fix the sign:
    %%%%%%%%%%%%%%%%%
    Eye = W_star/Wg;
    Find = find(sign(diag(Eye)) == -1);
    W_star(Find,:) = -W_star(Find,:);
    %% (7)Storing the frb estimates:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A_SymEfrb(:,:,b) = iD*W_star';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [gS, EdgS, dgS] = gdg(S,nonlin)

[d, n] = size(S);
gS = zeros(d,n);
dgS = zeros(d,n);
EdgS = zeros(d,1);

for i = 1:d
    g = nonlin;
    s = S(i,:);    
    if strcmp(g,'tanh')
        g = @(x) tanh(x);
        dg = @(x) 1-tanh(x).^2;
    elseif strcmp(g,'pow3')
        g = @(x) x.^3;
        dg = @(x) 3*x.^2;
    elseif strcmp(g,'gauss')
        g = @(x) x .* exp(- (1/2)*x.^2);
        dg =  @(x) (1 - x.^2) .* exp(- (1/2)*x.^2);
    elseif strcmp(g,'skew')
        g = @(x) x.^2;
        dg =  @(x) 0 ;
    elseif strcmp(g,'rt06')
        g = @(x) max(0,x-0.6).^2;
        dg = @(x) 2*max(0,x-0.6);
    elseif strcmp(g,'lt06')
        g = @(x) min(0,x+0.6).^2;
        dg = @(x) 2*min(0,x+0.6);
    elseif strcmp(g,'bt00')
        g = @(x) max(0,x).^2 - min(0,x).^2;
        dg = @(x) 2*max(0,x) - 2*min(0,x);
    elseif strcmp(g,'bt02')
        g = @(x) max(0,x-0.2).^2 - min(0,x+0.2).^2;
        dg = @(x) 2*max(0,x-0.2) - 2*min(0,x+0.2);
    elseif strcmp(g,'bt06')
        g = @(x) max(0,x-0.6).^2 - min(0,x+0.6).^2;
        dg = @(x) 2*max(0,x-0.6) - 2*min(0,x+0.6);
    elseif strcmp(g,'bt10')
        g = @(x) max(0,x-1.0).^2 - min(0,x+1.0).^2;
        dg = @(x) 2*max(0,x-1.0) - 2*min(0,x+1.0);
    elseif strcmp(g,'bt12')
        g = @(x) max(0,x-1.2).^2 - min(0,x+1.2).^2;
        dg = @(x) 2*max(0,x-1.2) - 2*min(0,x+1.2);
    elseif strcmp(g,'bt14')
        g = @(x) max(0,x-1.4).^2 - min(0,x+1.4).^2;
        dg = @(x) 2*max(0,x-1.4) - 2*min(0,x+1.4);
    elseif strcmp(g,'bt16')
        g = @(x) max(0,x-1.6).^2 - min(0,x+1.6).^2;
        dg = @(x) 2*max(0,x-1.6) - 2*min(0,x+1.6);
    elseif strcmp(g,'tan1')
        g = @(x) tanh(1.25*x);
        dg = @(x) 1.25*(1-tanh(1.25*x).^2);
    elseif strcmp(g,'tan2')
        g = @(x) tanh(1.5*x);
        dg = @(x) 1.5*(1-tanh(1.5*x).^2);
    elseif strcmp(g,'tan3')
        g = @(x) tanh(1.75*x);
        dg = @(x) 1.75*(1-tanh(1.75*x).^2);
    elseif strcmp(g,'tan4')
        g = @(x) tanh(2*x);
        dg = @(x) 2*(1-tanh(2*x).^2);
    elseif strcmp(g,'gau1')
        g = @(x) x .* exp(- (1.07/2)*x.^2);
        dg =  @(x) (1 - 1.07*x.^2) .* exp(- (1.07/2)*x.^2);
    elseif strcmp(g,'gau2')
        g = @(x) x .* exp(- (1.15/2)*x.^2);
        dg =  @(x) (1 - 1.15*x.^2) .* exp(- (1.15/2)*x.^2);
    elseif strcmp(g,'gau3')
        g = @(x) x .* exp(- (0.95/2)*x.^2);
        dg =  @(x) (1 - 0.95*x.^2) .* exp(- (0.95/2)*x.^2);
    else error('input g should be a string equal to tanh, pow3, gaus, lt06, rt06, bt00');
    end
    gS(i,:) = g(s);
    dgS(i,:) = dg(s);
    EdgS(i) = mean(dg(s));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[out]=bootrsp(in,B)
% out=bootrsp(in,B)
% Bootstrap resampling procedure.
% Inputs:
% in - input data
% B - number of bootstrap resamples (default B=1)
% Outputs:
% out - B bootstrap resamples of the input data
%
% For a vector input data of size [N,1], the resampling
% procedure produces a matrix of size [N,B] with columns
% being resamples of the input vector.
% For a matrix input data of size [N,M], the resampling
% procedure produces a 3D matrix of size [N,M,B] with
% out(:,:,i), i = 1,...,B, being a resample of the input
% matrix.
%
% Example:
%
% out=bootrsp(randn(10,1),10);
if (exist('B')~=1), B=1; end;
if (exist('in')~=1), error('Provide input data'); end;
s=size(in);
if length(s)>2,
error('Input data can be a vector or a 2D matrix only');
end;
if min(s)==1,
out=in(ceil(max(s)*rand(max(s),B)));
else
out=in(ceil(s(1)*s(2)*rand(s(1),s(2),B)));
end;