% hndWrtnS.m      
%            based on sftMx2LyrPt.m   23 Mar 2017
%             21 Apr 2017  
% Andrew P Paplinski, Monash University
% Hand-written digit recognition. MNIST dataset
% using stacked 28x28 images
% softmax 2-Layer NNet classifier with ReLU
% Pattern learning

clear
diary off
tic
init = 0; % 0 - continuation, 1 - initial run
czy = input('Initial run? y/n [y]: ', 's') ;
if (isempty(czy) || czy == 'y'), init = 1; end
if init == 1
    hdr = 'MNIST: hand-written digits classifier. Initial run.';
else
    hdr = 'MNIST: hand-written digits classifier. Continuation.';
end
% load MNIST/MNaniRT % X Y Xtst Ytst
load MNIST/mnistIL % Mimgs Mlbls MtstImgs MtstLbls
[r, c, N] = size(Mimgs) ;
X = reshape(Mimgs, r*c, N) ; % 784 by 60000
Y = Mlbls ;

% Limiting number of training patterns to:=========
N = 5000 ;
X = X(:, 1:N) ;
Y = Y(1:N) ;
% ================

[D, N] = size(X) ;     % dimensionality of X and the number of points
                       % D = 28, size of aniRt vectors                      
M = length(unique(Y)); % number of classes, dimensionality of the output Z
Xh = [X ; ones(1,N)] ; % adding bias input
% converting the expected class labels Y into outputs YY
YY = zeros(M,N) ;
for n = 1:N
    YY(Y(n)+1, n) = 1;  
end
 
% the result directory  
dr = 'rslts' ;   if isdir(dr) == 0,  mkdir(dr), end

%  setting up the time marker and the file_name to store the results 
tim = clock; dat = date ;
disp([hdr, date, sprintf('%4d:%02d',tim(4:5))])   
fnm = ['A',dat(4:6),sprintf('%02dh%02d%02d',tim(3:5))];
fid = fopen('CrntFnm.txt','wt');
fprintf(fid,'%s',fnm);  % writing a current base file name to CrntFnm.txt
fclose(fid);

dfn = [dr, '\', fnm] ; % the base file name for runs and maps
diary([dfn, '.txt']) 

fprintf(' Total number of training examples: %d', N)

% NNet specification:
% Dh hidden neurons with ReLU, M output neurons with Softmax
Dh = 40 ; % number of hidden neurons

lmda = 0; % 1e-1/N ; % lambda in reqularization
nepchs = 500  ;   % number of learning epochs
if init    
    % Weights initialization
    Wh = 0.01*(rand(Dh,D+1)-0.5) ; % 1st weight matrix is Dh by D+1
    Wz = 0.01*(rand(M,Dh+1)-0.5) ; % 2nd weight matrix is M by Dh+1
    % learning parameters
    etah = 1e-3 ; % 180/N ;    % learning rate (hidden)
    etaz = 1e-4 ; % 18/N ;     % learning rate (output)   
else 
    a = dir( 'rslts/*.mat') ;
    dfld = ['rslts/', a(end).name] ;
    load(dfld)   %  'Wh','Wz' , 'etah', 'etaz', 'L', 'LData'
    Lp = L; Lpdt = LData ;
end

fprintf('\n Number of learning epochs: %d', nepchs)
fprintf('\n etaz = %1.2e, etah = %1.2e', etaz, etah)

L = zeros(nepchs,1) ;     % total loss
LData = zeros(nepchs,1) ; % data loss
logpk = zeros (N,1) ;     % pattern log(pk_n)pk
ta = toc ;
fprintf('\n Time to complete in seconds:    ')

% learning
for epch=1:nepchs   % the epoch loop 
  for n=1:N % pattern loop, for each pattern  
    % evaluate all signals in the 2-layer Neural Network
    % first layer
    uh = Wh*Xh(:,n) ;    % postsynaptic activities in the hidden layer
    zh = [max(0,uh); 1]; % ReLU nonlinearity and bias
    % second layer
    % Calculating the postsynaptic activities (scores) M by 1
    uz = Wz*zh ; % each column (stimulus) gives the class scores  
                 % corresponding to the M classes 
    % the Softmax classifier -- the cross-entropy loss function
    ewx = exp(uz) ;     %  is M by 1
    s_ewx = sum(ewx) ;  %  is 1 by 1
    Pr = ewx./s_ewx ;   % probabilities, M by 1 = outputs z_n
           % selecting the correct probability for each class
    pk = sum(Pr.*YY(:,n)) ;  % is 1 by N , class probabilities
    logpk(n) = log(pk) ;
    
    % gradient of LData for the n-th pattern
    gL = Pr - YY(:,n) ;
    
    % The weight matrix update (including the regularization)
    dWz = gL*zh' + lmda*Wz ;  
     Wz = Wz - etaz*dWz ;

    % backpropagation from 2nd layer to the 1st layer
    delh = Wz(:,1:end-1)'*gL ;  % backprop error
    delh(uh<=0) = 0 ;           % backprop through ReLU
    % The weight matrix update (including the regularization)
    dWh = delh*Xh(:,n)' + lmda*Wh  ; 
    Wh = Wh - etah*dWh ;
  end   % pattern loop
  
  % data loss is an average cross entropy
  LData(epch) = -mean(logpk) ; % mean of -log(probs(pk))
  % lmda regularization loss
  Lreg = 0.5*lmda*(Wh(:)'*Wh(:) + Wz(:)'*Wz(:)) ;
  L(epch) = LData(epch) + Lreg  ;  % total loss
  if epch > 1 && L(epch) > L(epch-1) 
      fprintf('\n Rising L, epch = %d', epch)
      etaz = 0.75*etaz ;
      etah = 0.75*etah ;
      fprintf('\n etaz = %1.2e, etah = %1.2e', etaz, etah)
      fprintf('\n Time to complete in seconds:    ')
  end
  
  if rem(epch, 10) ==0
     % time to complete
     tmc = round((toc - ta)*(nepchs/epch - 1)) ;
     fprintf('\b\b\b\b %3d', tmc)
   end
    
end     %  the epoch loop  

fprintf('\n L(end) = %0.4f', L(end))
fprintf('\n Elapsed time is %d seconds\n', round(toc))

if init == 0
    L = [Lp; L] ;
    LData = [Lpdt; LData ] ;
end
% file name to store current network; from CrntFnm.txt
save(dfn,'Wh','Wz', 'etah', 'etaz', 'L', 'LData') %
fprintf('Saved as: %s.mat \n', dfn)

figure(1), clf
% plot(1:numel(L), [L LData], '-', numel(L), [L(end) LData(end)], 'o'), 
% legend('L', 'Ldata')
plot(1:numel(L), L, '-', numel(L), L(end), 'o'), 
% legend('L', 'Ldata')
grid on
title(['Loss function. Final value = ', num2str(L(end), '%0.4f')])
print('-f1', '-dpdf', dfn)

diary off

 