function [W,Cost,Grad] = ICA(X)

[N,T] = size(X); %X is the matrix of size NxT
[X, P] = pre_processing( X );

W = rand(N,N); %Initialize W
W = W/P;

max_iter = 500;
mu = 1 / 25;
min_cost = inf;
min_mu = 1/200;

kernels = zeros(N,max_iter);
Cost = zeros(max_iter,1);
time = zeros(max_iter,N);

for iter = 1 : max_iter
    entropy = zeros(N,1);
    Cost(iter) = -log(abs(det(W)));
    
    Wten = zeros(N-1,N,N);
    for kk = 1:N
        nout = [1:kk-1 kk+1:N]';
        Wten(:,:,kk) = W(nout, :);
    end
    
    for n = 1 : N
        
        Wn_last = Wten(:,:,n);
        
        temp1 = randn(N, 1);
        temp2 = Wn_last;
        h = temp1 - temp2'*((temp2*temp2')\temp2)*temp1;
        
        w = W(n,:)';
        y = w'*X;
        
        
        %%-----------------------------------------------------------------
        % Here we should define the entropy for each estimated source. This
        % estimation is based on the sei-cintinuous estimation technique.
        % Also we need to define phi which is the vector that contains the
        % derivatives. phi is a 1xT vector and entropy would be a single
        % number.
        %%-----------------------------------------------------------------
        
        grad = (X*phi')/T + h/(h'*w);
        
        
        grad = grad - (w'*grad)*w;
        
        
        grad_mon(n) = norm(grad);
        grad = grad / norm(grad);
        w1 = w + mu*grad;
        w1 = w1 / norm(w1);
        W_new(n,:) = w1';
        
        Grad(iter,n) = grad_mon(n);
        time(iter,n)=toc;
    end
    
    %termCriterion = max(1-abs(diag((W*P)*(W_new*P)')));
    oldwtchange = (W_new-W);
    delta=reshape(oldwtchange,1,N*N);
    termCriterion=delta*delta';
    termCriterion=norm(W_new-W,'fro');
    
    Grad(iter) = sum(grad_mon)/N^2;
    Cost(iter) = Cost(iter) + (N - sum(entropy));
    
   
    if Cost(iter) < min_cost
        cost_increase_counter = 0;
        min_cost = Cost(iter);
        W = W_new;
    else
        cost_increase_counter = cost_increase_counter + 1;
    end
        
    if cost_increase_counter > 10
        if mu > min_mu
            cost_increase_counter = 0;
            mu = mu/2;
            continue;
        else
            break;
        end
    else
        W = W_new;
    end
    
end


W = W*P;
end



function [X, P] = pre_processing( X )
% pre-processing program
[~,T] = size(X);
% remove DC
Xmean=mean(X,2);
X = X - Xmean*ones(1,T);
% spatio pre-whitening 1
R = X*X'/T;
P = inv_sqrtmH(R);  %P = inv(sqrtm(R));
X = P*X;
end


function A = inv_sqrtmH( B )
%
[V,D] = eig(B);
d = diag(D);
d = 1./sqrt(d);
A = V*diag(d)*V';
end
