% LEARNKURTOSIS	KURTOSIS learning rule.  
%
%  (c) Copyright by Ariel Tankus and Nathan Intrator  - Update Nov 19, 1995
%	
%		[dW,dB,risk] = LEARNKUTOSIS(ACTION,P,W,B,LP)
%		  ACTION - Either 'setup' to initialize the BCM parameters,
%		  	   or 'learn' to actually learn by KURTOSIS.
%		  P  - RxQ matrix of input vectors.
%		  W  - SxR matrix of weights.
%		  B  - Sx1 matrix of ouput vectors.
%		  LP,THETA,MOM4 - Learning Parameters. Affects only in 'setup'.
%		Returns:
%		  dW - a weight change matrix.
%		  dB - a bias change matrix.
%		  RISK - (optional) the risk (projection index).
%		  
%		[LP THETA MOM4] = LEARNKURTOSIS(ACTION)
%		  ACTION == 'getparams'. Get learning parameters.
%		Returns:
%		  LP - The learning parameters.
%		  
%		LP = [LR_BCM ETA]
%		  LR_BCM - Learning rate.
%		  ETA    - Feed Forward inhibition factor.
%		  THETA  - Dynamic threshold. If THETA == 0, it is
%			   computed based on the inputs.
%		  MOM4   - 4th momentum. Maintained for risk calculation.
%  Matlab version 5.3 changes, 2000 By Yuval Lifshitz.

%  Based on a BCM learning code.


function [dw,db,risk] = learnkurtosis(action,p,w,b,lp)

global k eta lr_kurtosis theta mom4 c sig_sup sig_inf

if (strcmp(action, 'getparams'))
   % dw is used to return parameters vector
   dw = [lr_kurtosis eta];
   db = theta;
   risk = mom4;
   return;
end


if (strcmp(action, 'setup'))
   % constants for BCM learning
   k       = 200; 			% memory length.
   lr_kurtosis  = lp(1);	% learning rate.
   eta     = lp(2); 		% inhibition factor.
   sig_sup = 5; 				% sigmoid supremum.
   sig_inf = -2; 			% sigmoid infimum.

end

% pack bias in weight and input matrices
w = [w b];
p = [p; ones(1, size(p, 2))];

c = w*p;
eta_vec = eta * ones(size(c, 1), 1);

% FF inhibited activity (called: ci - c inhibited)

N = (1 + eta) * c - eta_vec * sum(c);
ci = logsig(N);
ci = (sig_sup - sig_inf) * ci + sig_inf;

if (strcmp(action, 'setup'))
   if (isempty(theta))
      theta = mean(ci' .^ 2)';
   end
   if (isempty(mom4))

      mom4  = mean(ci' .^ 4)';
   end
   dw = 0;
   db = 0;
end


if (strcmp(action, 'learn'))
   theta = ((k-1) * theta + mean(ci' .^ 2)') ./ k;	% dynamic threshold.
   mom4 = ((k-1) * theta + mean(ci' .^ 4)') ./ k;		% 4th momentum.
   phi = 3.*ci .* (ci.^2 - theta*ones(1,size(ci,2)));		% modification function
   
   % derivative of sigmoid of ci
   dsig_ci = (sig_sup - sig_inf) * dlogsig(N,ci);
   
   % temporary results (before loop - optimization)
   tmp = phi .* dsig_ci;
   ones_row = ones(1,size(p,1));
   ones_col = ones(size(tmp,1),1);
   
   % Kurtosis_single: expectation of KURTOSIS non-inhibited output. 
   % i.e. treating each neuron as a single neuron.
   kurtosis_single = 0;
   for i = 1:size(p,2)
      kurtosis_single = kurtosis_single+(tmp(:,i) * ones_row) .* (ones_col * p(:,i)');
   end
   kurtosis_single = kurtosis_single ./ size(p,2);
   
   % take inhibition into account
   dw = lr_kurtosis * ( (1 + eta) * kurtosis_single - eta_vec * sum(kurtosis_single) );
   
   % extract bias
   db = dw(:, size(dw,2));
   dw = dw(:, 1:size(dw,2)-1);
end

% find the risk
if (nargout > 2)
   risk = -lr_kurtosis * sum(mom4 - 3.*(theta .^ 2));
end

