function w = ica(w,data) 
%w = ica(w,data)
%
% takes in an initial guess at the unmixing matrix (w) and a
% stream of data (mxn -- n signals of m samples)
% produces an umixing matrix
% (assumes all of the signals have zero mean and have been
% scaled to have a range of roughly 1)

[l, n] = size(data);
batchsize = floor(l/100);
eta = 0.5/batchsize; % eta needs to be inversely proportional
% to the batch size since we will be adding the sum of
% batchsize gradients all toether.
% note that eta must be much smaller for larger matrices
% and we must run the algorithm for longer
for s=1:20,
p = randperm(l);
data = data(p,:);
for tt=1:1000, % run through all of the data 10 times before
	% changing eta
	i = floor(rand(1)*(l-batchsize-1))+1;
	if (floor(tt/100)==tt/100),
		disp(s);
		disp(eta);
		disp(w);
	end;
	x = data(i:i+batchsize-1,:)';
	u = w*x;
	y = 1 ./ (1+exp(-u)); % compute non-linearity
	deltaw = (1 - 2*y)*u'*w; % second part of update rule
	% notice that the outer product automatically does the
	% sum over all of the examples for us
	deltaw = w*batchsize + deltaw;
	% the rule given was for one example (on-line learning)
	% So, since the term computed above was the sum over a 
	% batchsize number of examples, we also need the first
	% term to be the sum of batchsize number of examples.
	% Hence the w*batchsize instead of just w
	w = w + eta*deltaw;
end;
eta = eta*0.8; % this seems to work well for 2x2 matrices
end;
w = w'; % Bell & Sejnowski's algorithm assumes row vectors
% instead of column vectors
	
