function [W,Cost,mu,kernels,avetime,timematrix,Gradall] = ICA_EMK3_sid(X,ke)

[N,T] = size(X);
[X, P] = pre_processing( X );

% W = rand(N,N);
% W = W/P;
W = ICA_EBM(X); % Here we initialize using ICA-EBM. 

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

Params=struct('verbose',false);

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 % Here we can also use parfor. It works!
          
        Wn_last = Wten(:,:,n);
        
        temp1 = randn(N, 1);
        temp2 = Wn_last;
        h = temp1 - temp2'*((temp2*temp2')\temp2)*temp1;
        
        w = W(n,:)';
        y = w'*X;
        
        
        [~, Lambda, paramLocalConst, ~]=EMK_mod2(y,'numLocalConst',ke,'greedySearch',0); %This 3 corresponds to the number of local kernels 
        lparam = paramLocalConst{1,1};
        len = length(Lambda)-4;
        kernels(n,iter) = (length(Lambda) - 4);
        
        ALPHA2 = [];
        if len > 0
            MU = repmat(lparam(1:len,1), 1, T);
            SIGMA = repmat(lparam(1:len,2), 1, T);
            Y = repmat(y, len, 1);
            ALPHA2 = sum(kern_val(Y, MU, SIGMA), 2) / T;
        end
        yty = y*y';
        ALPHA = [1; mean(y); yty/T; mean(y./(1+y.^2)); ALPHA2];
        
        entropy(n) = Lambda' * ALPHA;
        
        
        phi = zeros(1,T);
        for i = 1:length(Lambda)
            if i <=4
                phi = phi + Lambda(i) * rg_dif(y,i);
            else
                phi = phi + Lambda(i) * kern_dif(y, lparam(i-4,1), lparam(i-4,2));
            end
        end
        
        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';
        
        Gradall(iter,n) = grad_mon(n);
        %time(iter,n)=toc;
    end 
    
    %Grad(iter) = sum(grad_mon)/N^2; Maybe here we need to take the
    %weighted average in order to take into account the number of kernels
    Grad(iter) = norm(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
    
    if Params.verbose
        fprintf('\n Step %d:, Cost: %f, mu: %f, Grad: %f', iter,Cost(iter), mu, Grad(iter));
    end % options.verbose
    
end

W = W*P;
timeserial = sum(time,2);
avetime = mean(timeserial);
timematrix = time;
end

function val = rg_dif(y, method)
switch method
    case 1
        val = zeros(size(y));
    case 2
        val = ones(size(y));
    case 3
        val = 2.*y;
    case 4
        val = (1-y.^2)./((1+y.^2).^2);
end
end

function val = kern_val(y, mu, sigma)
val = exp(-((y - mu).^2)./(2*(sigma).^2));
end

function val = kern_dif(y, mu, sigma)
val = -(y - mu)./(sigma^2).*exp(-((y - mu).^2)./(2*sigma^2));
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
