function [w] = trace;

%
% implementation of Foldiak's (1991) original trace rule paper for Bi/CNS
% homework set 221
%
% return weight matrix after simulation
%
% Wolfgang Einhaeuser-Treyer
% May, 2006
%

% graphical output on/off
verboseInput  = 0;  % switch on if you want to see the input in figure 1 (slow)
verboseMatrix = 1;  % switch off to make simulation faster

% plot matix every ... sweeps
plotEvery = 10;

% parameters
delta    = 0.2;  % trace parameter (see original paper)
alpha    = 0.02; % learning rate
gridSize = 8;    % size of SC grid
nCC      = 4;    % number of complex cells
nOri     = 4;    % number of orientations / this is fixed

nSweeps     = 500;

% simple cells 
x = zeros(nOri,gridSize,gridSize);

% complex cells
y = zeros(1,nCC);

% trace
y_trace = zeros(1,nCC);

% weights
w = 0.1*rand(nCC,nOri,gridSize,gridSize); % post / pre

% main simulation
for sweepC = 1 : nSweeps
  
  %
  % generate input for current sweep
  %

  fprintf('\r%4.4d',sweepC)
  
  % pick orientation and direction for current sweep
  sweepOri  = ceil(nOri*rand(1));
  direction = (rand(1)>0.5);
  
  %
  % diagonal is longer
  %
  
  if (direction)
	if sweepOri>2
	  tVec = [2*gridSize-1:-1:1];
	else
	  tVec = [gridSize:-1:1];
	end  
  else
	if sweepOri>2
	  tVec = [1:2*gridSize-1];
	else
	  tVec = [1:gridSize];
	end  
  end
  
  
  %
  % perform the sweep
  %
  
  for t = tVec
	
	% generate simple cell activity
	x = zeros(nOri,gridSize,gridSize);

	switch sweepOri
	  
	case 1 
	  % horizontal
	  x(sweepOri,:,t) = 1;
	case 2
	  % vertical
	  x(sweepOri,t,:) = 1;
	case 3
	  % slanted
	  if t <= gridSize
		onY = gridSize-[1:t]+1;
		onX = [t:-1:1];
	  else
		
		tNew = t-gridSize+1;
		
		onY = gridSize-[tNew:gridSize]+1;
		onX = [gridSize:-1:tNew];
	  end
	  onI = sub2ind([gridSize,gridSize],onY,onX);
	  x(sweepOri,onI) = 1;

	case 4
	    % back slanted
	  if t <= gridSize
		onX = [1:t];
		onY = [t:-1:1];
	  else
		
		tNew = t-gridSize+1;
		onX = [tNew:gridSize];
		onY = [gridSize:-1:tNew];
	  end
	  onI = sub2ind([gridSize,gridSize],onY,onX);
	  x(sweepOri,onI) = 1;

	end

	
	%
	% show input
	%
	
	if verboseInput
	  % plot simple cell activity
	  figure(1);
	  for l = 1 : nOri
		subplot(1,nOri,l);
		imagesc(squeeze(x(l,:,:)));
		switch l	
		case 1 
		  title('SC horizontal');
		case 2
		  title('SC vetical');
		case 3
		  title('SC slanted');
		case 4
		  title('SC back-slanted');
		end
		drawnow;
	  end
	end

	
	
	%
	% compute output (postsynaptic/CC) activity
	%

	% linear sum
	for i = 1 : nCC
	  y(i) = w(i,:)*x(:);
	end
	
	% apply maximum rule
	y = double(y==max(y));
	
	
	%
	% LEARNING (this can be done more effectively in matlab without
	% loops of course, but for clarity I stick to the notation in 
	% the paper
	%

	for i = 1 : nCC
	  
	  % collapse the presynaptic dimensions from 3 to 1 (j)
	  % fortunately, matlab does this correctly

	  for j = 1 : nOri*gridSize*gridSize
		w(i,j) = ...
			w(i,j) + ...
			alpha*y_trace(i)*(x(j)-w(i,j));
	  end
	end
	
	
	
	% update trace
	y_trace = (1-delta)*y_trace+delta*y;
  
  
	
  end
  
  
  if (verboseMatrix)&(mod(sweepC,plotEvery)==0)
	figure(2);
	clf;
	for i = 1 : nCC
	  
	  for j1 = 1 : nOri
		
		subplot(nCC,nOri,nOri*(i-1)+j1);
		imagesc(squeeze(w(i,j1,:,:)),[min(w(:)),max(w(:))]);
		title(sprintf('CC %d / SC ori: %d',i,j1));
	  end
	end
	
	drawnow;
  end
  
end