function [patterns, targets, W] = ICA(patterns, targets, params)

%Reshape the data points using the independent component analysis algorithm
%Inputs:
%	train_patterns	- Input patterns
%	train_targets	- Input targets
%	params			- [Output dimension, Learning rate]
%
%Outputs
%	patterns		- New patterns
%	targets			- New targets
%	W				- Reshape martix

[r,c]			= size(patterns);
[dimension, eta] = process_params(params);

if (r < dimension),
   error('Output dimension cannot be larger than the input dimension')
end

%Whiten the data to zero mean and unit covariance
patterns    = patterns - mean(patterns')'*ones(1,c);
[v, d]	    = eig(cov(patterns',1));
Aw			= v*inv(sqrtm(d));
patterns    = Aw'*patterns;

%Move data to the range of [-1,1]
patterns = (patterns - min(patterns')'*ones(1,c))./((max(patterns') - min(patterns'))'*ones(1,c));
patterns = patterns*2-1;

%Find the weight matrix
W			= randn(r);
for i = 1:c,
   y		= patterns(:,i);
   phi	    = activation(y);
   dW		= eta*(eye(r) - phi*y')*W;
   
   %Break if algorithm diverges
   if (max(max(dW)) > 1e3),
       disp(['Algorithm diverged after ' num2str(i) ' iterations'])
       break
   end
   
   W		= W + dW;   
   
   %If the algorithm converged, exit
   if (max(max(W)) < 1e-2),
       disp(['Algorithm converged after ' num2str(i) ' iterations'])
       break
   end
end

%Take only the most influential outputs
power		= sum(abs(W)');
[m, in]	= sort(power);
W			= W(in(r-dimension+1:r),:);

%Calculate new patterns
patterns = W*patterns;


%End ICA

function phi = activation(y)
%Activation function for ICA
%phi=(3/4)*y.^11+(25/4)*y.^9+(-14/3)*y.^7+(-47/4)*y.^5+(29/4)*y.^3;
phi = 0.5*y.^5 + 2/3*y.^7 + 15/2*y.^9 + 2/15*y.^11 - 112/3*y.^13 + 128*y.^15 - 512/3*y.^17;