function [A, W] = debugica(X, whiteningMatrix, dewhiteningMatrix, approach, ...
			numOfIC, g, finetune, a1, a2, myy, stabilization, ...
			epsilon, maxNumIterations, maxFinetune, initState, ...
			guess, sampleSize, displayMode, displayInterval, ...
			s_verbose);
% DebugICA - Debug main algorithm of Fixed point ICA.
% Perform independent component analysis using Hyvarinen's fixed point
% algorithm. Outputs an estimate of the mixing matrix A and its inverse W.
%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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

b_verbose=1;
approachMode = 2;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

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

switch lower(stabilization)
 case 'on'
  stabilizationEnabled = 1;
 case 'off'
  if myy ~= 1
    stabilizationEnabled = 1;
  else
    stabilizationEnabled = 0;
  end
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;
  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
  end
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFLATION APPROACH
if approachMode == 2
  B = zeros(vectorSize);
  tt = 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;
    
        % Take a random initial vector of lenght 1 and orthogonalize it
        % with respect to the other vectors.
        if initialStateMode == 0
            w = rand(vectorSize, 1) - .5;
        elseif initialStateMode == 1
            w=whiteningMatrix*guess(:,round);
        end
        w = w - B * B' * 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
              figure(3); 
              tt(:, round) = w;
              temp = X'*tt;
              icaplot('dispsig',temp(:,1:numOfIC)');
              drawnow;
              figure(round+3);
              hist(temp(:,round)');
              figure(3);
              pause;  
            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
            end
%                                     B(:, round) = w;
%                                     temp = X'*B;
%                                  icaplot('dispsig',temp(:,1:numOfIC)');
%                                  drawnow;
%                                    pause; 
            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 = (Xub * ((Xub' * 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
      
      % 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)
  A = real(A);
  W = real(W);
end

