% TBPKURTOSIS2 Train 2-layer feed-forward network w/backpropagation & KURTOSIS.
%
%  (c) Copyright by Ariel Tankus and Nathan Intrator  - Update Nov 19, 1995
%
% [W1,B1,W2,B2,TE,TR,RISK_REC] = TBP2(W1,B1,F1,W2,B2,F2,P,T,TP, ...
%	,TEST_P,TEST_T)
%   Wi - SixR weight matrix of ith layer.
%   Bi - Six1 bias vector of ith layer.
%   Fi - Transfer function (string) of ith layer.
%   P  - RxQ matrix of input vectors.
%   T  - S2xQ matrix of target vectors.
%   TP - Training parameters (optional).
%   THETA - Init value for dynamic threshold for kurtosis (0 = calc from input).
%   MOM4 - Init value for 4th momentum for kurtosis (0 = calc from input).
%   TEST_P - Input test sample set.
%   TEST_T - Target test sample set.
% Returns:
%   Wi - new weights.
%   Bi - new biases.
%   TE - the actual number of epochs trained.
%   TR - training record: [row of errors for bp]
%   RISK_REC - risk record: [row of errors for kurtosis]
%
% Training parameters are:
%   TP(1)  - Epochs between updating display, default = 25.
%   TP(2)  - Epoch number to start training, default = 1.
%   TP(3)  - Epoch number to stop training, default = 1000.
%   TP(4)  - Sum-squared error goal, default = 0.02.
%   TP(5)  - Learning rate for BP, default = 0.01.
%   TP(6)  - Learning rate for Kurtosis, default = 0.01.
%   TP(7)  - Eta, the inhibition factor for Kurtosis, default = 0.02.
%   TP(8)  - the weight of BP relative to Kurtosis (range: [0,1]), default = 0.9.
%   TP(9), TP(10) - The dimensions of input image, default = 1.
%   TP(11) - Number of successive epochs w/o classification error.
%		If stable<=0, don't stop training by CE, default = 1.
%
% Missing parameters and NaN's are replaced with defaults.
% 
%	    
%	Mark Beale, 1-31-92
%	Revised 12-15-93, M.B.
%	Copyright (c) 1992-94 by the MathWorks, Inc.
%	$Revision: 1.1 $  $Date: 1994/01/11 16:29:32 $
%
%	Modified by Ariel Tankus, 1995.

%  Matlab version 5.3 changes, 2000 By Yuval Lifshitz.

%  Based on a BCM training code.

function [w1,b1,w2,b2,i,tr,risk_rec]=tbpkurtosis2(w1,b1,f1,w2,b2,f2,p,t,tp, ...
	test_P,test_T)


if (nargin < 8)

   error('Not enough arguments.');

end

% should test set be checked for error rate?
test_flag = (nargin >= 12);

% TRAINING PARAMETERS
if (nargin == 7)

   tp = [];

end


tp = nndef(tp,[25 1 1000 0.02 0.01 0.01 0.02 0.9 1000 [1, size(w1,2)] 1]);
df     = tp(1);
se     = tp(2);
me     = tp(3);
eg     = tp(4);
lr_bp  = tp(5);
lr_kurtosis = tp(6);
eta    = tp(7);
bpw    = tp(8); 					% BP Weight.
dim    = [tp(9) tp(10)]; 		% image/weights dimensions.
stable = tp(11);
df1    = feval(f1,'deriv');
df2    = feval(f2,'deriv');

% PRESENTATION PHASE

z1  = w1*p;

q   = size(z1,2);
a1  = feval(f1,netsum(z1, concur(b1,q)));

z2  = w2*a1;

q   = size(z2,2);
a2  = feval(f2,netsum(z2, concur(b2,q)));
e   = t-a2;
MSE = mean(mean(e.^2));
[CLERR_PERCENT clerr] = classerr(t, a2);

% INIT KURTOSIS
lp = [lr_kurtosis eta];
[dw1_kurtosis, db1_kurtosis, risk] = learnkurtosis('setup', p, w1, b1, lp);

% TRAINING RECORD
tr 	= zeros(1,me);
tr(1) = MSE;
risk_rec 	= zeros(1,me);
risk_rec(1) = risk;

% Test Set Error Rate
if (test_flag)

   test_z1  = w1*test_P;

   test_q	  = size(test_z1,2);
   test_a1  = feval(f1,netsum(test_z1, concur(b1,test_q)));

   test_z2  = w2*test_a1;

   test_q	  = size(test_z2,2);
   test_a2  = feval(f2,netsum(test_z2, concur(b2,test_q)));
   test_err = classerr(test_T, test_a2);
end  

% CE - Classification Error rate (on training set)
if (test_flag)
   message = sprintf('BPKURTOSIS: %%g/%g, MSE=%%g, RISK=%%g, dW=%%g, CE=%%g%%%%, TEST=%%g%%%%, kurtosis/bp=%%.3g\n',me);
else
   message = sprintf('BPKURTOSIS: %%g/%g epochs, MSE=%%g, RISK=%%g, dW=%%g, CE=%%g%%%%, kurtosis/bp=%%.3g\n',me);
end

fprintf(message,se - 1,MSE,risk,0,CLERR_PERCENT,test_err,0)

zero_err_counter = 0;

for i=se : me
   
   % CHECK PHASE
   if (CLERR_PERCENT == 0 & stable > 0)
      % `zero_err_counter' counts the number of successive epochs with
      % 0% misclassification rate.
      zero_err_counter = zero_err_counter + 1;
      
      % stop training if 0% error for `stable' epochs
      if (zero_err_counter >= stable)
         break;
      end

      
   else
      zero_err_counter = 0;
   end
   
   if (MSE < eg)

      i=i-1;

      break;

   end
   
   % BACKPROPAGATION PHASE
   d2 = feval(df2,e,a2);

   
   d1 = feval(df1,d2,a1);

   
   % LEARNING PHASE

   
   % BP learning
   if (bpw ~= 0)
      [dw1_bp,db1_bp] = learnbp(p,d1,lr_bp);
      [dw2,db2] = learnbp(a1,d2,lr_bp);
   end
   
   % KURTOSIS learning
   if (bpw ~= 1)
      [dw1_kurtosis,db1_kurtosis,risk] = learnkurtosis('learn',p,w1,b1);
   end
   
   % Combining the two learning methods

   if (bpw == 1)
      dw1 = dw1_bp;
      db1 = db1_bp;
   elseif bpw == 0
      dw1 = dw1_kurtosis;
      db1 = db1_kurtosis;
   else
      dw1 = bpw * dw1_bp + (1-bpw) * dw1_kurtosis;
      db1 = bpw * db1_bp + (1-bpw) * db1_kurtosis;
   end
   
   w1 = w1 + dw1; b1 = b1 + db1;
   w2 = w2 + dw2; b2 = b2 + db2;
   
   % PRESENTATION PHASE

   z1  = w1*p;

   q   = size(z1,2);
   a1  = feval(f1,netsum(z1, concur(b1,q)));

   z2  = w2*a1;

   q   = size(z2,2);
   a2  = feval(f2,netsum(z2, concur(b2,q)));
   e   = t-a2;
   MSE = mean(mean(e.^2));
   [CLERR_PERCENT clerr] = classerr(t, a2);
   
   % TRAINING RECORD
   tr(i+1) = MSE;
   risk_rec(i+1) = risk;
   
   if (rem(i,df) == 0)
      
      % Test Set Error Rate
      if (test_flag)

         test_z1  = w1*test_P;

         test_q	  = size(test_z1,2);
         test_a1  = feval(f1,netsum(test_z1, concur(b1,test_q)));

         test_z2  = w2*test_a1;

         test_q	  = size(test_z2,2);
         test_a2  = feval(f2,netsum(test_z2, concur(b2,test_q)));
         test_err = classerr(test_T, test_a2);
      end	
      % PLOTTING
      bp_norm = norm(bpw*dw1_bp);

      if (bp_norm < 1e-5)

         kurtosisbp_ratio = inf;

      else

         kurtosisbp_ratio = norm((1 - bpw)*dw1_kurtosis) / bp_norm;

      end

      
      fprintf(message,i,MSE,risk,norm(dw1_kurtosis)/lr_kurtosis,CLERR_PERCENT, ...
         test_err,kurtosisbp_ratio)
      
   end
   

end

% TRAINING RECORD
tr = tr(1:(i+1));
risk_rec = risk_rec(1:(i+1));

if (rem(i,df) ~= 0)
   % Test Set Error Rate
   if (test_flag)

      test_z1  = w1*test_P;

      test_q	  = size(test_z1,2);
      test_a1  = feval(f1,netsum(test_z1, concur(b1,test_q)));

      test_z2  = w2*test_a1;

      test_q	  = size(test_z2,2);
      test_a2  = feval(f2,netsum(test_z2, concur(b2,test_q)));
      test_err = classerr(test_T, test_a2);
   end
   
   % PLOTTING

   bp_norm = norm(bpw*dw1_bp);

   if bp_norm < 1e-5

      bp_norm = 1e-5 ;

   end
   kurtosisbp_ratio = norm((1-bpw)*dw1_kurtosis) / bp_norm;
   fprintf(message,i,MSE,risk,norm(dw1_kurtosis)/lr_kurtosis,CLERR_PERCENT, ...
      test_err,kurtosisbp_ratio)

end

% WARNINGS
if (MSE > eg)
   disp(' ')
   disp('TRAINBPKURTOSIS: Network error did not reach the error goal.')
   disp('  Further training may be necessary, or try different')
   disp('  initial weights and biases and/or more hidden neurons.')
   disp(' ')
end

