function [A, W] = fpica(X, whiteningMatrix, dewhiteningMatrix, approach, ...
    numOfIC, g, finetune, a1, a2, myy, stabilization, ...
    epsilon, maxNumIterations, maxFinetune, initState, ...
    guess, sampleSize, displayMode, displayInterval, ...
    s_verbose);
%FPICA - Fixed point ICA. Main algorithm of FASTICA.
%
% [A, W] = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach,
%        numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon,
%        maxNumIterations, maxFinetune, initState, guess, sampleSize,
%        displayMode, displayInterval, verbose);
%
% Perform independent component analysis using Hyvarinen's fixed point
% algorithm. Outputs an estimate of the mixing matrix A and its inverse W.
%
% whitesig                              :the whitened data as row vectors
% whiteningMatrix                       :can be obtained with function whitenv
% dewhiteningMatrix                     :can be obtained with function whitenv
% approach      [ 'symm' | 'defl' ]     :the approach used (deflation or symmetric)
% numOfIC       [ 0 - Dim of whitesig ] :number of independent components estimated
% g             [ 'pow3' | 'tanh' |     :the nonlinearity used
%                 'gaus' | 'skew' ]
% finetune      [same as g + 'off']     :the nonlinearity used in finetuning.
% a1                                    :parameter for tuning 'tanh'
% a2                                    :parameter for tuning 'gaus'
% mu                                    :step size in stabilized algorithm
% stabilization [ 'on' | 'off' ]        :if mu < 1 then automatically on
% epsilon                               :stopping criterion
% maxNumIterations                      :maximum number of iterations
% maxFinetune                           :maximum number of iteretions for finetuning
% initState     [ 'rand' | 'guess' ]    :initial guess or random initial state. See below
% guess                                 :initial guess for A. Ignored if initState = 'rand'
% sampleSize    [ 0 - 1 ]               :percentage of the samples used in one iteration
% displayMode   [ 'signals' | 'basis' | :plot running estimate
%                 'filters' | 'off' ]
% displayInterval                       :number of iterations we take between plots
% verbose       [ 'on' | 'off' ]        :report progress in text format
%
% EXAMPLE
%       [E, D] = pcamat(vectors);
%       [nv, wm, dwm] = whitenv(vectors, E, D);
%       [A, W] = fpica(nv, wm, dwm);
%
%
% This function is needed by FASTICA and FASTICAG
%
%   See also FASTICA, FASTICAG, WHITENV, PCAMAT

% @(#)$Id: fpica.m,v 1.7 2005/06/16 12:52:55 jarmo Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variable for stopping the ICA calculations from the GUI
global g_FastICA_interrupt;
if isempty(g_FastICA_interrupt)
    clear global g_FastICA_interrupt;
    interruptible = 0;
else
    interruptible = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values

if nargin < 3, error('Not enough arguments!'); end
[vectorSize, numSamples] = size(X);
if nargin < 20, s_verbose = 'on'; end
if nargin < 19, displayInterval = 1; end
if nargin < 18, displayMode = 'on'; end
if nargin < 17, sampleSize = 1; end
if nargin < 16, guess = 1; end
if nargin < 15, initState = 'rand'; end
if nargin < 14, maxFinetune = 100; end
if nargin < 13, maxNumIterations = 1000; end
if nargin < 12, epsilon = 0.0001; end
if nargin < 11, stabilization = 'on'; end
if nargin < 10, myy = 1; end
if nargin < 9, a2 = 1; end
if nargin < 8, a1 = 1; end
if nargin < 7, finetune = 'off'; end
if nargin < 6, g = 'pow3'; end
if nargin < 5, numOfIC = vectorSize; end     % vectorSize = Dim
if nargin < 4, approach = 'defl'; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the data

if ~isreal(X)
    error('Input has an imaginary part.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for verbose

switch lower(s_verbose)
    case 'on'
        b_verbose = 1;
    case 'off'
        b_verbose = 0;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for approach

switch lower(approach)
    case 'symm'
        approachMode = 1;
    case 'defl'
        approachMode = 2;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''approach''\n', approach));
end
if b_verbose, fprintf('Used approach [ %s ].\n', approach); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for numOfIC

if vectorSize < numOfIC
    error('Must have numOfIC <= Dimension!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the sampleSize
if sampleSize > 1
    sampleSize = 1;
    if b_verbose
        fprintf('Warning: Setting ''sampleSize'' to 1.\n');
    end
elseif sampleSize < 1
    if (sampleSize * numSamples) < 1000
        sampleSize = min(1000/numSamples, 1);
        if b_verbose
            fprintf('Warning: Setting ''sampleSize'' to %0.3f (%d samples).\n', ...
                sampleSize, floor(sampleSize * numSamples));
        end
    end
end
if b_verbose
    if  b_verbose & (sampleSize < 1)
        fprintf('Using about %0.0f%% of the samples in random order in every step.\n',sampleSize*100);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for nonlinearity.

switch lower(g)
    case 'pow3'
        gOrig = 10;
    case 'tanh'
        gOrig = 20;
    case {'gaus', 'gauss'}
        gOrig = 30;
    case 'skew'
        gOrig = 40;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''g''\n', g));
end
if sampleSize ~= 1
    gOrig = gOrig + 2;
end
if myy ~= 1
    gOrig = gOrig + 1;
end

if b_verbose,
    fprintf('Used nonlinearity [ %s ].\n', g);
end

finetuningEnabled = 1;
switch lower(finetune)
    case 'pow3'
        gFine = 10 + 1;
    case 'tanh'
        gFine = 20 + 1;
    case {'gaus', 'gauss'}
        gFine = 30 + 1;
    case 'skew'
        gFine = 40 + 1;
    case 'off'
        if myy ~= 1
            gFine = gOrig;
        else
            gFine = gOrig + 1;
        end
        finetuningEnabled = 0;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''finetune''\n', ...
            finetune));
end

if b_verbose & finetuningEnabled
    fprintf('Finetuning enabled (nonlinearity: [ %s ]).\n', finetune);
end

switch lower(stabilization)
    case 'on'
        stabilizationEnabled = 1;
    case 'off'
        if myy ~= 1
            stabilizationEnabled = 1;
        else
            stabilizationEnabled = 0;
        end
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''stabilization''\n', ...
            stabilization));
end

if b_verbose & stabilizationEnabled
    fprintf('Using stabilized algorithm.\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some other parameters
myyOrig = myy;
% When we start fine-tuning we'll set myy = myyK * myy
myyK = 0.01;
% How many times do we try for convergence until we give up.
failureLimit = 5;


usedNlinearity = gOrig;
stroke = 0;
notFine = 1;
long = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for initial state.

switch lower(initState)
    case 'rand'
        initialStateMode = 0;
    case 'guess'
        if size(guess,1) ~= size(whiteningMatrix,2)
            initialStateMode = 0;
            if b_verbose
                fprintf('Warning: size of initial guess is incorrect. Using random initial guess.\n');
            end
        else
            initialStateMode = 1;
            if size(guess,2) < numOfIC
                if b_verbose
                    fprintf('Warning: initial guess only for first %d components. Using random initial guess for others.\n', size(guess,2));
                end
                guess(:, size(guess, 2) + 1:numOfIC) = ...
                    rand(vectorSize,numOfIC-size(guess,2))-.5;
            elseif size(guess,2)>numOfIC
                guess=guess(:,1:numOfIC);
                fprintf('Warning: Initial guess too large. The excess column are dropped.\n');
            end
            if b_verbose, fprintf('Using initial guess.\n'); end
        end
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''initState''\n', initState));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for display mode.

switch lower(displayMode)
    case {'off', 'none'}
        usedDisplay = 0;
    case {'on', 'signals'}
        usedDisplay = 1;
        if (b_verbose & (numSamples > 10000))
            fprintf('Warning: Data vectors are very long. Plotting may take long time.\n');
        end
        if (b_verbose & (numOfIC > 25))
            fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
        end
    case 'basis'
        usedDisplay = 2;
        if (b_verbose & (numOfIC > 25))
            fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
        end
    case 'filters'
        usedDisplay = 3;
        if (b_verbose & (vectorSize > 25))
            fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
        end
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''displayMode''\n', displayMode));
end

% The displayInterval can't be less than 1...
if displayInterval < 1
    displayInterval = 1;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if b_verbose, fprintf('Starting ICA calculation...\n'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SYMMETRIC APPROACH
if approachMode == 1,
    
    % set some parameters more...
    usedNlinearity = gOrig;
    stroke = 0;
    notFine = 1;
    long = 0;
    
    A = zeros(vectorSize, numOfIC);  % Dewhitened basis vectors.
    if initialStateMode == 0
        % Take random orthonormal initial vectors.
        B = orth (randn (vectorSize, numOfIC));
    elseif initialStateMode == 1
        % Use the given initial vector as the initial state
        B = whiteningMatrix * guess;
    end
    
    BOld = zeros(size(B));
    BOld2 = zeros(size(B));
    
    % This is the actual fixed-point iteration loop.
    for round = 1:maxNumIterations + 1,
        if round == maxNumIterations + 1,
            fprintf('No convergence after %d steps\n', maxNumIterations);
            fprintf('Note that the plots are probably wrong.\n');
            if ~isempty(B)
                % Symmetric orthogonalization.
                B = B * real(inv(B' * B)^(1/2));
                
                W = B' * whiteningMatrix;
                A = dewhiteningMatrix * B;
            else
                W = [];
                A = [];
            end
            return;
        end
        
        if (interruptible & g_FastICA_interrupt)
            if b_verbose
                fprintf('\n\nCalculation interrupted by the user\n');
            end
            if ~isempty(B)
                W = B' * whiteningMatrix;
                A = dewhiteningMatrix * B;
            else
                W = [];
                A = [];
            end
            return;
        end
        
        
        % Symmetric orthogonalization.
        B = B * real(inv(B' * B)^(1/2));
        
        % Test for termination condition. Note that we consider opposite
        % directions here as well.
        minAbsCos = min(abs(diag(B' * BOld)));
        minAbsCos2 = min(abs(diag(B' * BOld2)));
        
        if (1 - minAbsCos < epsilon)
            if finetuningEnabled & notFine
                if b_verbose, fprintf('Initial convergence, fine-tuning: \n'); end;
                notFine = 0;
                usedNlinearity = gFine;
                myy = myyK * myyOrig;
                BOld = zeros(size(B));
                BOld2 = zeros(size(B));
                
            else
                if b_verbose, fprintf('Convergence after %d steps\n', round); end
                
                % Calculate the de-whitened vectors.
                A = dewhiteningMatrix * B;
                break;
            end
        elseif stabilizationEnabled
            if (~stroke) & (1 - minAbsCos2 < epsilon)
                if b_verbose, fprintf('Stroke!\n'); end;
                stroke = myy;
                myy = .5*myy;
                if mod(usedNlinearity,2) == 0
                    usedNlinearity = usedNlinearity + 1;
                end
            elseif stroke
                myy = stroke;
                stroke = 0;
                if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
                    usedNlinearity = usedNlinearity - 1;
                end
            elseif (~long) & (round>maxNumIterations/2)
                if b_verbose, fprintf('Taking long (reducing step size)\n'); end;
                long = 1;
                myy = .5*myy;
                if mod(usedNlinearity,2) == 0
                    usedNlinearity = usedNlinearity + 1;
                end
            end
        end
        
        BOld2 = BOld;
        BOld = B;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Show the progress...
        if b_verbose
            if round == 1
                fprintf('Step no. %d\n', round);
            else
                fprintf('Step no. %d, change in value of estimate: %.3g \n', round, 1 - minAbsCos);
            end
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Also plot the current state...
        switch usedDisplay
            case 1
                if rem(round, displayInterval) == 0,
                    % There was and may still be other displaymodes...
                    % 1D signals
                    icaplot('dispsig',(X'*B)');
                    drawnow;
                end
            case 2
                if rem(round, displayInterval) == 0,
                    % ... and now there are :-)
                    % 1D basis
                    A = dewhiteningMatrix * B;
                    icaplot('dispsig',A');
                    drawnow;
                end
            case 3
                if rem(round, displayInterval) == 0,
                    % ... and now there are :-)
                    % 1D filters
                    W = B' * whiteningMatrix;
                    icaplot('dispsig',W);
                    drawnow;
                end
            otherwise
        end
        
        switch usedNlinearity
            % pow3
            case 10
                B = (X * (( X' * B) .^ 3)) / numSamples - 3 * B;
            case 11
                % optimoitu - epsilonin kokoisia eroja
                % t�m� on optimoitu koodi, katso vanha koodi esim.
                % aikaisemmista versioista kuten 2.0 beta3
                Y = X' * B;
                Gpow3 = Y .^ 3;
                Beta = sum(Y .* Gpow3);
                D = diag(1 ./ (Beta - 3 * numSamples));
                B = B + myy * B * (Y' * Gpow3 - diag(Beta)) * D;
            case 12
                Xsub=X(:, getSamples(numSamples, sampleSize));
                B = (Xsub * (( Xsub' * B) .^ 3)) / size(Xsub,2) - 3 * B;
            case 13
                % Optimoitu
                Ysub=X(:, getSamples(numSamples, sampleSize))' * B;
                Gpow3 = Ysub .^ 3;
                Beta = sum(Ysub .* Gpow3);
                D = diag(1 ./ (Beta - 3 * size(Ysub', 2)));
                B = B + myy * B * (Ysub' * Gpow3 - diag(Beta)) * D;
                
                % tanh
            case 20
                hypTan = tanh(a1 * X' * B);
                B = X * hypTan / numSamples - ...
                    ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / numSamples * ...
                    a1;
            case 21
                % optimoitu - epsilonin kokoisia
                Y = X' * B;
                hypTan = tanh(a1 * Y);
                Beta = sum(Y .* hypTan);
                D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
                B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
            case 22
                Xsub=X(:, getSamples(numSamples, sampleSize));
                hypTan = tanh(a1 * Xsub' * B);
                B = Xsub * hypTan / size(Xsub, 2) - ...
                    ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / size(Xsub, 2) * a1;
            case 23
                % Optimoitu
                Y = X(:, getSamples(numSamples, sampleSize))' * B;
                hypTan = tanh(a1 * Y);
                Beta = sum(Y .* hypTan);
                D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
                B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
                
                % gauss
            case 30
                U = X' * B;
                Usquared=U .^ 2;
                ex = exp(-a2 * Usquared / 2);
                gauss =  U .* ex;
                dGauss = (1 - a2 * Usquared) .*ex;
                B = X * gauss / numSamples - ...
                    ones(size(B,1),1) * sum(dGauss)...
                    .* B / numSamples ;
            case 31
                % optimoitu
                Y = X' * B;
                ex = exp(-a2 * (Y .^ 2) / 2);
                gauss = Y .* ex;
                Beta = sum(Y .* gauss);
                D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
                B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
            case 32
                Xsub=X(:, getSamples(numSamples, sampleSize));
                U = Xsub' * B;
                Usquared=U .^ 2;
                ex = exp(-a2 * Usquared / 2);
                gauss =  U .* ex;
                dGauss = (1 - a2 * Usquared) .*ex;
                B = Xsub * gauss / size(Xsub,2) - ...
                    ones(size(B,1),1) * sum(dGauss)...
                    .* B / size(Xsub,2) ;
            case 33
                % Optimoitu
                Y = X(:, getSamples(numSamples, sampleSize))' * B;
                ex = exp(-a2 * (Y .^ 2) / 2);
                gauss = Y .* ex;
                Beta = sum(Y .* gauss);
                D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
                B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
                
                % skew
            case 40
                B = (X * ((X' * B) .^ 2)) / numSamples;
            case 41
                % Optimoitu
                Y = X' * B;
                Gskew = Y .^ 2;
                Beta = sum(Y .* Gskew);
                D = diag(1 ./ (Beta));
                B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
            case 42
                Xsub=X(:, getSamples(numSamples, sampleSize));
                B = (Xsub * ((Xsub' * B) .^ 2)) / size(Xsub,2);
            case 43
                % Uusi optimoitu
                Y = X(:, getSamples(numSamples, sampleSize))' * B;
                Gskew = Y .^ 2;
                Beta = sum(Y .* Gskew);
                D = diag(1 ./ (Beta));
                B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
                
            otherwise
                error('Code for desired nonlinearity not found!');
        end
    end
    
    
    % Calculate ICA filters.
    W = B' * whiteningMatrix;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Also plot the last one...
    switch usedDisplay
        case 1
            % There was and may still be other displaymodes...
            % 1D signals
            icaplot('dispsig',(X'*B)');
            drawnow;
        case 2
            % ... and now there are :-)
            % 1D basis
            icaplot('dispsig',A');
            drawnow;
        case 3
            % ... and now there are :-)
            % 1D filters
            icaplot('dispsig',W);
            drawnow;
        otherwise
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFLATION APPROACH
if approachMode == 2
    
    B = zeros(vectorSize);
    
    % The search for a basis vector is repeated numOfIC times.
    round = 1;
    
    numFailures = 0;
    
    while round <= numOfIC,
        myy = myyOrig;
        usedNlinearity = gOrig;
        stroke = 0;
        notFine = 1;
        long = 0;
        endFinetuning = 0;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Show the progress...
        if b_verbose, fprintf('IC %d ', round); end
        
        % Take a random initial vector of lenght 1 and orthogonalize it
        % with respect to the other vectors.
        if initialStateMode == 0
            w = randn (vectorSize, 1);
        elseif initialStateMode == 1
            w=whiteningMatrix*guess(:,round);
        end
        w = w - B * B' * w; %orthogonalization - B is the previous w
        w = w / norm(w);
        
        wOld = zeros(size(w));
        wOld2 = zeros(size(w));
        
        % This is the actual fixed-point iteration loop.
        %    for i = 1 : maxNumIterations + 1
        i = 1;
        gabba = 1;
        while i <= maxNumIterations + gabba
            if (usedDisplay > 0)
                drawnow;
            end
            if (interruptible & g_FastICA_interrupt)
                if b_verbose
                    fprintf('\n\nCalculation interrupted by the user\n');
                end
                return;
            end
            
            % Project the vector into the space orthogonal to the space
            % spanned by the earlier found basis vectors. Note that we can do
            % the projection with matrix B, since the zero entries do not
            % contribute to the projection.
            w = w - B * B' * w;
            w = w / norm(w);
            
            if notFine
                if i == maxNumIterations + 1
                    if b_verbose
                        fprintf('\nComponent number %d did not converge in %d iterations.\n', round, maxNumIterations);
                    end
                    round = round - 1;
                    numFailures = numFailures + 1;
                    if numFailures > failureLimit
                        if b_verbose
                            fprintf('Too many failures to converge (%d). Giving up.\n', numFailures);
                        end
                        if round == 0
                            A=[];
                            W=[];
                        end
                        return;
                    end
                    % numFailures > failurelimit
                    break;
                end
                % i == maxNumIterations + 1
            else
                % if notFine
                if i >= endFinetuning
                    wOld = w; % So the algorithm will stop on the next test...
                end
            end
            % if notFine
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Show the progress...
            if b_verbose, fprintf('.'); end;
            
            
            % Test for termination condition. Note that the algorithm has
            % converged if the direction of w and wOld is the same, this
            % is why we test the two cases.
            if norm(w - wOld) < epsilon | norm(w + wOld) < epsilon
                if finetuningEnabled & notFine
                    if b_verbose, fprintf('Initial convergence, fine-tuning: '); end;
                    notFine = 0;
                    gabba = maxFinetune;
                    wOld = zeros(size(w));
                    wOld2 = zeros(size(w));
                    usedNlinearity = gFine;
                    myy = myyK * myyOrig;
                    
                    endFinetuning = maxFinetune + i;
                    
                else
                    numFailures = 0;
                    % Save the vector
                    B(:, round) = w;
                    
                    % Calculate the de-whitened vector.
                    A(:,round) = dewhiteningMatrix * w;
                    % Calculate ICA filter.
                    W(round,:) = w' * whiteningMatrix;
                    
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Show the progress...
                    if b_verbose, fprintf('computed ( %d steps ) \n', i); end
                    
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Also plot the current state...
                    switch usedDisplay
                        case 1
                            if rem(round, displayInterval) == 0,
                                % There was and may still be other displaymodes...
                                % 1D signals
                                temp = X'*B;
                                icaplot('dispsig',temp(:,1:numOfIC)');
                                drawnow;
                            end
                        case 2
                            if rem(round, displayInterval) == 0,
                                % ... and now there are :-)
                                % 1D basis
                                icaplot('dispsig',A');
                                drawnow;
                            end
                        case 3
                            if rem(round, displayInterval) == 0,
                                % ... and now there are :-)
                                % 1D filters
                                icaplot('dispsig',W);
                                drawnow;
                            end
                    end
                    % switch usedDisplay
                    break; % IC ready - next...
                end
                %if finetuningEnabled & notFine
            elseif stabilizationEnabled
                if (~stroke) & (norm(w - wOld2) < epsilon | norm(w + wOld2) < ...
                        epsilon)
                    stroke = myy;
                    if b_verbose, fprintf('Stroke!'); end;
                    myy = .5*myy;
                    if mod(usedNlinearity,2) == 0
                        usedNlinearity = usedNlinearity + 1;
                    end
                elseif stroke
                    myy = stroke;
                    stroke = 0;
                    if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
                        usedNlinearity = usedNlinearity - 1;
                    end
                elseif (notFine) & (~long) & (i > maxNumIterations / 2)
                    if b_verbose, fprintf('Taking long (reducing step size) '); end;
                    long = 1;
                    myy = .5*myy;
                    if mod(usedNlinearity,2) == 0
                        usedNlinearity = usedNlinearity + 1;
                    end
                end
            end
            
            wOld2 = wOld;
            wOld = w;
            
            switch usedNlinearity
                % pow3
                case 10
                    w = (X * ((X' * w) .^ 3)) / numSamples - 3 * w;
                case 11
                    EXGpow3 = (X * ((X' * w) .^ 3)) / numSamples;
                    Beta = w' * EXGpow3;
                    w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
                case 12
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    w = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2) - 3 * w;
                case 13
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    EXGpow3 = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2);
                    Beta = w' * EXGpow3;
                    w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
                    % tanh
                case 20
                    hypTan = tanh(a1 * X' * w);
                    w = (X * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / numSamples;
                case 21
                    hypTan = tanh(a1 * X' * w);
                    Beta = w' * X * hypTan;
                    w = w - myy * ((X * hypTan - Beta * w) / ...
                        (a1 * sum((1-hypTan .^2)') - Beta));
                case 22
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    hypTan = tanh(a1 * Xsub' * w);
                    w = (Xsub * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / size(Xsub, 2);
                case 23
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    hypTan = tanh(a1 * Xsub' * w);
                    Beta = w' * Xsub * hypTan;
                    w = w - myy * ((Xsub * hypTan - Beta * w) / ...
                        (a1 * sum((1-hypTan .^2)') - Beta));
                    % gauss
                case 30
                    % This has been split for performance reasons.
                    u = X' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    w = (X * gauss - sum(dGauss)' * w) / numSamples;
                case 31
                    u = X' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    Beta = w' * X * gauss;
                    w = w - myy * ((X * gauss - Beta * w) / ...
                        (sum(dGauss)' - Beta));
                case 32
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    u = Xsub' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    w = (Xsub * gauss - sum(dGauss)' * w) / size(Xsub, 2);
                case 33
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    u = Xsub' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    Beta = w' * Xsub * gauss;
                    w = w - myy * ((Xsub * gauss - Beta * w) / ...
                        (sum(dGauss)' - Beta));
                    % skew
                case 40
                    w = (X * ((X' * w) .^ 2)) / numSamples;
                case 41
                    EXGskew = (X * ((X' * w) .^ 2)) / numSamples;
                    Beta = w' * EXGskew;
                    w = w - myy * (EXGskew - Beta*w)/(-Beta);
                case 42
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    w = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
                case 43
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    EXGskew = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
                    Beta = w' * EXGskew;
                    w = w - myy * (EXGskew - Beta*w)/(-Beta);
                    
                otherwise
                    error('Code for desired nonlinearity not found!');
            end
            
            
% TEST % % % % % % % % % % % % % % % % % % % % % % % % %
interval = 1;
if mod(i,interval)==0
    modifyA =1; %for now...
    nx = 20;
    ny = 20;
    sigfix = 3.7/4; % sig of PSF...
    
    if exist ('modifyA','var')
        if modifyA == 1
            %                     fprintf ('Modifying matrix conevrgence of w...')
            Atemp = dewhiteningMatrix * w;
%             signum = 1;
%             if abs(min(Atemp))>max(Atemp) %or sum(Atemp)<0...?
%                 signum = -1;
%             end
%             if abs(min(Atemp))>max(Atemp) %or sum(Atemp)<0...?
%                 Atemp = -Atemp;
%             end
            absA = abs(Atemp);
            maxA = max(absA);
            
            gf = (gauss2d([nx ny],[nx/2,ny/2],sigfix,maxA));
            Aconv = (conv2(gf, reshape (Atemp, nx, ny),'same')).^5;
            Atemp = reshape(Aconv,ny*ny,1);
            
%             [x_mu, y_mu, sig] = fitgauss2d(reshape(absA, nx, ny),sigfix); %fit to abs...
%             Atemp = (Atemp + signum * reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1))/2;
            % Atemp = reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1);
        end
    end
    w = whiteningMatrix*Atemp;
end
            
% TEST % % % % % % % % % % % % % % % % % % % % % % % % %

% %%% TEST version old %%%%%%%%%
% if mod(i,5)==0
%     modifyA =1; %for now...
%     nx = 20;
%     ny = 20;
%     sigfix = 3.7/4; % sig of PSF...
%     
%     if exist ('modifyA','var')
%         if modifyA == 1
%             %                     fprintf ('Modifying matrix conevrgence of w...')
%             Atemp = dewhiteningMatrix * w;
%             if abs(min(Atemp))>max(Atemp) %or sum(Atemp)<0...?
%                 Atemp = -Atemp;
%             end
%             [x_mu, y_mu, sig] = fitgauss2d(reshape(Atemp, nx, ny),sigfix);
%             maxA = max(abs(Atemp));
%             Atemp = (Atemp + reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1))/2;
%             % Atemp = reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1);
%         end
%     end
%     w = whiteningMatrix*Atemp;
% 
% %%% TEST version old %%%%%%%%%

% Normalize the new w.
            w = w / norm(w);    
            i = i + 1;
        end
        round = round + 1;
    end
    if b_verbose, fprintf('Done.\n'); end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Also plot the ones that may not have been plotted.
    if (usedDisplay > 0) & (rem(round-1, displayInterval) ~= 0)
        switch usedDisplay
            case 1
                % There was and may still be other displaymodes...
                % 1D signals
                temp = X'*B;
                icaplot('dispsig',temp(:,1:numOfIC)');
                drawnow;
            case 2
                % ... and now there are :-)
                % 1D basis
                icaplot('dispsig',A');
                drawnow;
            case 3
                % ... and now there are :-)
                % 1D filters
                icaplot('dispsig',W);
                drawnow;
            otherwise
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In the end let's check the data for some security
if ~isreal(A)
    if b_verbose, fprintf('Warning: removing the imaginary part from the result.\n'); end
    A = real(A);
    W = real(W);
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunction
% Calculates tanh simplier and faster than Matlab tanh.
function y=tanh(x)
y = 1 - 2 ./ (exp(2 * x) + 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Samples = getSamples(max, percentage)
Samples = find(rand(1, max) < percentage);
